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Abstract 

By the use of elegant operational properties of Legendre 
wavelets, a direct computational algorithm for evaluation 
the optimal control and trajectory of linear singular 
systems with a quadratic cost functional is developed. 
The state variable, state- rate, and the control vector are 
expanded in Legendre orthonormal basis with 
unknown coefficient. The technique is straightforward 
and very convenient for computer programming. An 
illustrative example is given to demonstrate the 
applicability of the proposed method. 

Keywords: Optimal control; singular systems; Legendre 
wavelets; operational matrix; orthogonal functions 

1. Introduction 

Wavelet theory is relatively new and is an emerging area 
in mathematical research. In recent years, wavelets have 
found their way into different fields of science and 
engineering. Wavelets permit the accurate representation 
of a variety of functions and operators. Orthogonal 
functions and polynomial series have been received 
considerable attention in dealing with various problems 
of dynamical systems. The main characteristic of this 
technique is that it reduces these problems to those of 
solving a system of algebraic equations, thus greatly 
simplifying the problems. The approach in [1] is based on 
converting the underlying differential equation into 
integral equations through integration, approximating 
various signals involved in the equation by orthogonal 
series and using the operational matrix of integration, to 
eliminate the integral operations. Special attention has 
been given to application of the Legendre wavelets [2], 
the linear Legendre wavelets [1], and the cosine and sine 
(CAS) wavelets which are used to solve Lredholm 
integral equation in [3]. Recently, constrained optimal 
control problems had been discussed many authors 
[4, 5, 6], the advantage of using Legendre wavelets as 
compared to the other wavelets is that in the present 
method the operational matrix are zeros, and integration 



of the product of two Legendre wavelet function vectors 
is an identity matrix, hence making Legendre wavelets 
computationally very attractive. Singular systems, also 
commonly called generalized or descriptor systems in the 
literature appear in many practical situations including 
engineering systems, economic systems, network 
analysis, and biological systems (see e.g. [7, 8, 9]). In 
fact, many systems in the real life are singular in nature. 
They are usually simplified as or approximated by non- 
singular models because there is still lacking of efficient 
tools to tackle problems related to such systems. The 
structural analysis of linear singular systems, using either 
algebraic or geometric approach, has attracted 
considerable attention from many researchers during the 
last three decades (see e.g. [10, 11, 12, 13, 14, 15, 16, and 
17] and the references cited therein). Generally speaking, 
almost all the research works dealing with singular 
systems are the natural extension of those results for non- 
singular counterparts, although it is much harder in 
obtaining solutions associated with singular systems. In 
this paper we extend the Legendre wavelets technique to 
study the problem of optimal control of linear singular 
systems with quadratic performance. The state 
variable x(^) , state rate x(^)and control variable 
u{t ) are expanded by Legendre wavelets with unknown 
coefficients. The advantage of this method is that we can 
choose series numbers willfully. 

Also the properties of Legendre wavelets are used to 
relate the unknown coefficients of x{t) to the 

coefficients of x(t) .Using this method the dynamics of 

the systems are converted to a system of algebraic 
equations .the necessary conditions of optimality are 
imposed to solve the quadratic programming by using the 
Lagrange multipliers. 

The reminder of the paper is organized as follows: 
section (2) is about the preliminaries and problem 
statement; section (3) focuses on the main result, section 
(4) demonstrate a benchmark example. Finally section (5) 
is the conclusion. 
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2. Preliminaries and Problem Statement 

2.1 Properties of the Legendre wavelets 

2.1.1 Multiresolution of analysis (MRA) 

In recent years, wavelets have found their way into many 
different fields of science and engineering. Wavelets 
constitute a family of functions constructed from dilation 
and translation of single function called the mother 
wavelet. When the dilation parameter a and the 
translation parameter b vary continuously, we have the 
following family of continuous wavelets [18]. 

i i— t — b 

y/ ab (0 = H 2 )’ a,b e R, 0 

a 

If we restrict the parameters a and b to discrete values 

as a = 2~ k ,b = n2~ k , then 

k 

Vkn(. t ) = 2 2 W(? k t-n) 



form an orthogonal basis. MRA has been used to 
construct the wavelets. We need the following definition: 

the increasing sequence {Vj} eZ of subset of L 2 (R) 

with scaling function (j) is called MRA if it satisfies the 
condition 

(/) 'Uj Vj Is dense in L 2 (R) , 

(«) n/ f { 0}, 

iiii) f(t)sVjOf(2- J t)sV 0 , 

(/v) {(p{t — n)} nGZ Is an orthogonal basis forP^ . 



Using the definition of MRA, the quence 

l 

{2 2 <f>(2 J t-n)} neZ forms an orthonormal basis for 

v j ■ 

Let (//(t)be the mother wavelet, then we have [19] 
y/(t) = '^a(n)</>(2t - n) 

n&Z 

and by suitable conditions, the sequence 

k 

{2 2 \ff(2 k t — n)} k n&z forms an orthonormal basis 
for L 2 ( R ) . 

Legendre wavelets y/ m n (t) = i// (k ,n,m ,t) have 
four arguments; k = 1,2,3, . . . , 



n = 2n -1 , n = 1,2,3, . . .,2* l 9 m is the order for 
Legendre polynomials and t is the normalized time. 
They are defined on interval [0, 1) by: 



( 1 ) 



( m + L m (2 k t-h), ^ 
0 



1 fi + 1 
2 k 

oterwise 



Where m = 0,1,2, ...,M — 1 ,n = 1,2,3, ... ,2^ ! . In Eq. 

1 - 

(1), the coefficient ( m H — ) 2 is for orthonormality. 



Here, L m ( t ) are the well-known Legendre polynomials 
of order m which are orthogonal with respect to the 
weight function w(/) = 1 and satisfy the following 
recursive formula: 

Z 0 (7) = 1,A (t) = t, 



L m+l (t) = 



(2m + \\ T r 

77 “ tL m (0 ~ 

V m + 1 



m 



m + 1 



4»-i(0>w = 1,2,3,... 



2.1.2. Function approximation 

Suppose a signal f(t) has finite energy over the 
interval [0,1) . Then f(t) can be represented by Legendre 
wavelets as follows 



OO 00 

= ( 2 ) 

n=\ m = 0 

where c n,m in which (, ) denotes the 

inner product. 

If the infinite series in Eq. (2) is truncated, then Eq. (2) 
can be written as 



/(o*ZZ c «,^(o=c T m (3) 

n - 1 m - 0 

Where C and 4 > (?) are 2 k Mxl matrices given by 
C = [c w ,c n ,.. ,C W _J,C 20 ,. . ;C 2M _ J,. . .,C 2 t_, 0 ,. . ;C 2 k-i M _ j] 



T (4) 

, y/ 2M _, ( 0 , • • • , ( 0 , • • • , ¥ 2 ^ m _ x ( 0 ] 

The integration of the vector ) defined in Eq. (4) 
can be obtained as 

l'¥(t')dt' = P'¥(t), (5) 

Where P is the (2 k ~ l M) x (2 k ~ l M) operational 
matrix for integration as 




L 


H 


H 


... H 


0 


L 


H 


... H 




0 












H 


0 


0 




0 L 



( 6 ) 
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H = 



And 



2 0 
0 0 

0 0 






1 



A 

3 

0 



^00 

a/3 

0 4 0 
3/5 

A o A 

5/3 5/7 

0 o 4 0 

7a/ 5 



0 0 0 0 



0 0 0 0 



0 



J2M-A 



(2M-3){2M-l 

0 



The integration of product of two Legendre wavelets 
vector function is obtained as: 



d=i = 

Jo 

Where / is an identity matrix. 



(7) 



x N (t) = C T N '¥(t) 



(13) 



Where 

A n — [^M0 5 a NU 9 • * * 9 a N\M -\ 9 ^ N 20 9 * * * 9 a N2M -\ 9 * * * 

’ ^ N 2 k ~ l 0 9 * * * 9 a N2 k - x M-\ ] 

(14) 

tfz = 10 9 ^Zll’***’^Z 1M-1 9 ^Z 20 9 * ’ * 9 ^ Z2M-X 9 * * * 

, U Z2 k-i 0 , . . . , ^ Z 2 k ~ l M-\ ] 

(15) 

^ V — 10 9 1 5 • • • 9 C N 1M-1 9 C N 20 9 * * * 9 C N 2M-1 9 * * * 

? 2 k ~ l 0 9 * * * 9 C N 2 k ~ l M-\ ] 

(16) 

Where {cL Nij } , {u Zij . } and {c Nij . } are unknown. 

Let 

(Ax(t)) N =Y^(t) (17) 

(Bu(t)) N = 4^(0 (18) 

(19) 



(Ex(t)) N =W£'¥(t) 

Using equation (17)-(19) for each 
N,N = 1,2,39 ^ the equation (8) has the form 

ff;^)=(Y;u;mo 



2.2 Problem statement 

Consider linear time invariant (LTI) singular system 

Ex{t) = v4x(£) + 5w(£) (8) 

x(0)=x 0 (9) 



Where x e R p is the state of the system, u <= R 2 is the 

input vector of the system, E G R pxp is assumed to be 
singular with 0 < rank ( E ) < p ,and A,B are real 
constant matrices of appropriate size. The problem is to 
find the optimal control u(t ) and the corresponding state 

trajectory x(t\ 0 < t < Z, satisfying equation (8) and 
(9) while minimizing the cost function 

./ = A(L)SxiL) + 1 £ [x J (t)Qx(t) + u T (t)Ru(t) ]dt (10) 

In equation (10), T denotes transposition, S,Q and R 
are matrices of appropriate dimensions, S and Q are 
symmetric positive semi-defmite matrices and R is a 
symmetric positive definite matrix. 



3. Main result 

3.1. Approximation of the singular system 

By expanding X N {t\ U z (t) and x N ( t ) in Legendre 

wavelet functions, we determine the following 
approximated values, i.e. for N = 1,2,3, . . . n, and 



Z = l,2,3,...z we get 

% (o = 4^(0 



u z {t) = U T z ^{t) 



( 11 ) 

( 12 ) 



By equating the coefficients of the same- order Legendre 
wavelets, we obtain for N = 1,2,3, . . .n and 

W T N ={Yl+L\) (20) 

Also, we can have 

x(t ) = J x(t')dt' + x(0) 

4^(0 = c ^(0 + a :^(0 

Equivalently 

A t n -C t n P-K t n = 0 (21) 

K n = [x(0),0, . . . ,0|x(0),0, . . . ,0|. . .|x(0),0, . . . ,0] T 

We define 

2 k ~ l M 

F \s ~ a NS ~ ^j C Nr^rS ~ ^ NS (^2) 

r = 1 

N = l,2,...,n , S = 1,2,...,2 W M 

3.2. The performance index approximation 

We now approximate the performance index J by using 
Legendre wavelets. 

Let 






10 ^11 ••• 1 ^20 ••• ••• a \i~\) ■“ a \t A M 



V \ 



(23) 




3 



Automatic Control and System Engineering Journal, ISSN 1687-4811, Volume 9, Issue II, ICGST LLC, Delaware, USA, December 2009 





(^10 


^ 1 1 • • • H IM-l 


U\ 20 • • • 


^YM-X • • • 


“irt - J T 


p= 


( 








\t 




j^ZIO 


U Z\ 1 • • • U Z\M-\ 


U Z20 • • 


• U Z2M-l • 


•• V-'o ••• u z^-'m-\) J 






(¥\t) 


> 






{ m= 






II 

<9? 








K 


v J (o / 







(24) 



(25) 



Note that a,/? and (t) md'¥ 2 (t) are matrices of 

order N2 M x 1, Z2 k ~ l M x 1 and N x 7V2^ M and 

Z x Z2 A 1 M respectively. Using equation (11) and (12) 
and (23) - (25) the state and control vector can be 
expressed as 

x(t) = % (t)a (26) 

u(t) = W 2 (t)j3 (27) 

Substituting equation (26) and (27) in equation (10) we 
get 

j =ia T 'P 1 T (l)S'P 1 (l)« 







a 



Y 2 T (t)R¥ 2 (t)dt 



P 



(28) 



Equation (28) can be computed more efficiently by using 
equation (7) and writhing as 

J =U t ['P(1)'P t (1)®5]«+U t (Z) ®Q)a (29) 

+^P\D®R)P 

Where ® denotes the kronecker product [20] 

Remark: As it is shown in (29) the dynamic cost 
function is converted to a statistic cost function and it 
simplifies the problem considerably. 



dL 

92 



(32) 



Remark: As it shown in (6) the operational matrix P 
has so many zeros and this fact results in a considerably 
saving of computing time. 



4. Illustrative example 



Consider the linear singular system 


[21] 


1 


"o r 


ii 


= 


"1 o~ 


*1 


+ 


'O' 


0 0 _ 


x 2 




0 1 


x 2 




1 



X(0) = 




(33) 

(34) 



With the cost function 

J = l£ 2 [x T (r)X(r) + C/ 2 (r) ]dr ( 35 ) 

The problem is to find the optimal control u(t) which 

minimizes equation (35) subject to (33) and (34). 

To obtain the Legendre wavelets approximations, we first 
covert the interval [0,2] to [0,1] by means of the 



transformation^ = — . We determine Legendre wavelets 

approximations for M = 5, k = 2, 

Let 

A = [hio a ill a U2 a m a U4 a i20 a \2l a i22 a \23 ^124 ] 

A * [^210 ^211 a 2\2 a 2U a 2\4 a 220 a 22\ a 222 a 223 a 224 ] 

u = \u 10 U n U l2 w l3 U u u 2 o U 21 u 2 2 U 2 3 W 24 ] 

Q = [ C 110 C lll C 112 C 113 C 114 C 120 C 121 C 122 C 123 C 124 ] 

f*2 = [ C 210 C 211 C 212 C 213 C 214 C 220 C 221 C 222 C 223 C 224 ] 

Applying equations (28), we obtain the following 
approximation for J 



min 



j=Y2L<m+ a L 

n = 1 m=0 



+ U 



2 

nm 



3.3. Solution of the optimization problem 

The optimal control problem has been reduced to a 
parameter optimization problem which can be stated as 
follows. Lind a and /? so that J(a,j3) is minimized 
subject to the constraints in equations (20) - (22) 

Let 

n 2 k ~ x M 

L = J + X ( 3 °) 

N = 1 S = 1 

Where A, ns* ^ = 1,2 5 = 1,2,. . .,2* _1 Af are 

unknown Lagrange multipliers. Then the necessary 
conditions for extremum are 



dL 

dd NS 

dL 

dU 7 , 



= 0 



N = 1,2 Z = 1,2,... , z , 



(31) 



= 0 



S = 1,2,... ,2 k ~ l M : 



Using equations (30) and (31) the Legendre wavelets 
approximation for X x {t) and u(t ) can be calculated. 
In table I and TI , a comparison is made between the 
exact solution and approximated values of Xj(^)and 
u(t ) using present method for M = 5,k = 2, and the 
exact solution. 

5. Conclusion 

The Legendre wavelets and the associated matrix of 
integration (operational matrix) are applied to solve the 
optimal control of singular systems. The method is based 
upon reducing a linear singular quadratic optimization 
problem to a set of algebraic equations. The integration 
of the product of two Legendre wavelet function vectors 
is an identity matrix, hence making Legendre wavelets 
computationally very attractive. Illustrative example is 
included to demonstrate the validity and applicability of 
the technique. 
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table I . Estimated and exact value of x (t) 



X, (0 


Solution 

number 


Time 


Exact 

solution 


Legendre 

wavelets 


Error 


1 


0.0 


1.0000 


1.0012 


0.0012 


2 


0.1 


0.7536 


0.7536 


0.0000 


3 


0.2 


0.5680 


0.5671 


0.0009 


4 


0.3 


0.4280 


0.4273 


0.0007 


5 


0.4 


0.3226 


0.3225 


0.0001 


6 


0.5 


0.2431 


0.2446 


0.0015 


7 


0.6 


0.1832 


0.1832 


0.0000 


8 


0.7 


0.1381 


0.1379 


0.0002 


9 


0.8 


0.1041 


0.1039 


0.0002 


10 


0.9 


0.0784 


0.0784 


0.0000 


11 


1.0 


0.0591 


0.0595 


0.0004 



table n . Estimated and exact value of u(t) 
u(t ) 



Solution 

number 


Time 


Exact 

solution 


Legendre 

wavelets 


Error 


1 


0.0 


0.7071 


0.7079 


0.0008 


2 


0.1 


0.5329 


0.5328 


0.0001 


3 


0.2 


0.4016 


0.4010 


0.0006 


4 


0.3 


0.3027 


0.3021 


0.0006 


5 


0.4 


0.2281 


0.2280 


0.0001 


6 


0.5 


0.1719 


0.1729 


0.0010 


7 


0.6 


0.1296 


0.1296 


0.0000 


8 


0.7 


0.0976 


0.0975 


0.0001 


9 


0.8 


0.0736 


0.0735 


0.0001 


10 


0.9 


0.0555 


0.0554 


0.0001 


11 


1.0 


0.0418 


0.0420 


0.0002 
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Abstract 

The major incidents are rare but their potential consequences 
can be very serious. The regional or national black out in the 
system causes the paralysis of the economic activity and an 
upheaval of the social life all the more marked as the 
breakdown is longer and wider. It is necessary for this reason 
to underline the correlation between the geographical extent 
of a black out zone and the duration of its realimentation, 
being given the specific difficulties of the reconstitution of a 
great system after collapse. In the case this occur, a plan of 
restoration of the network must be established beforehand. It 
is thus justified to make the provisions suitable to control the 
incidents when they happen, in order to limit the 
consequences in terms of extent, duration of the black out 
and to prepare a fast and effective realimentation for the 
affected consumers. 

The approach of the majority of the electricity companies on 
these problems resulted in putting a whole of measurements 
called plan of defence, which based on the electrotechnical 
phenomena brought during incidents and their own dynamics. 
In this paper, we present the principle of defence plan design 
into preventive and curative action. Implementation of the 
defence plans integrating two groups of defensive 
measurements : local measurements of defence 

recommended generally for the lines of interconnection and 
centralized measurements of defence or co-ordinated called 
special protection system (SPS). A synthesis study of 1 1 1 
centralized defence Plans made it possible to raise the strong 
and weak points of the whole of these defence plans. Finally 
the defence plans of the maghreb electrical system , as well 
as the encountered difficulties of implementation, were 
presented, commented and evaluated in order to create a new 
future optimal strategy. 

Key words: collapse of the electrical system, strategy of defence, 
curative defence , preventive defence. 

1. Introduction 

The analysis of the continuity of service introduced the 
enumeration of standard states characterizing the operation 
of the high voltage electrical system. The disturbances 
affecting the system can result in more degraded states. The 
actions of defence make it possible to bring back the system 
to less degraded states. 



To restore the state of the system primarily means to preserve 
the integrity of the network. This is to say: preserve as far as 
possible the maximum of network equipment and the 
interconnection grid system in service, to bring back the 
basic variables, the voltage, the frequency, the load of the 
connections in their respective limits, in order to avoid 
additional releases; but also to preserve and maintain 
available the manufacturing units, connected or not to the 
system. This restoration has a price: the shedding of load, the 
insolation of zones of the initial incident or victims of 
consecutive disturbances. 

The defence measures ensure on the one hand the detection 
of abnormal situations and on the other hand the setting of 
automatic actions, generally preestablished, allowing to 
preserve the integrity of the network, under acceptable 
conditions of performance. 

In the current state of the technique, and in accordance with 
the tradition in protection or adjustment, the establishment 
rests preferentially on local measurements and criteria, at 
least for the ultimate criterion of release of the 
countermeasure process. These conditions are indeed 
regarded the base of the robustness. A tendency today is to 
promote the use of systematic information coming from 
remote areas benefitting from the channels of information 
transfer high flow which furrow the networks. 

For each mechanism, it is possible to define actions of 
defence, and for each one of them to draw up a list of the 
means likely to stop the development of the degradation 
process of the system’s state. These means must be used in a 
radical way according to the process in which they fit. The 
design of the protection plans of such a network depends on 
the structure of the network topology, the site of 
production... etc. The paper is organized as follows; the 1st 
section after this introduction includes the design of the 
defense plans, the 2nd sections deals with the 

implementation of the defense plan. Section 3 includes the 
protection plan in the Maghreb electrical system. Section 4 
states the improvement activities and the future 

recommendations. And finally the conclusion of the paper. 
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2. Design of defence in depth 

At a major or generalized incident, various physical 
phenomena of degradation appear and follow one another. 
The safety of the system rests on the setting to be put a 
provisions of various nature, adapting to the dynamics of 
each phenomenon and preventing, detecting and treating the 
disfunctions being able to lead to their emergence and/or to 
control their evolution. These provisions are called lines of 
defence. The setting of the successive lines of defence 
constitutes the concept of defence in-depth (preventive and 
curative). 

When a disturbance affects an element of network, the 
elementary system of protection starts, and if the network 
were in N-l safety, its return towards a stable operation point 
is normally assured. So that an initially banal incident turns 
into a major incident, a series of dysfunctions must intervene 
such as a bad adjustment of protections, a refusal of opening 
breaker or that the point of initial operation does not 
guarantee a sufficient operation margin, etc. 

The extension of the incident can also occur by various 
mechanisms directly affecting the state, the appearance of 
static instabilities, imbalance between production and 
consumption resulting from a parcelling out of the network, 
instabilities of voltage resulting from the absence of reserves 
or the automatic operation of tap changer of the distribution 
transformers, the inadequate load shedding, or the reaction 
speed of the operators unsuited to incident dynamics. 

In any general information, the disturbances (causes) bring 
back the electric system from an initial state (F 0 ) to another 
state characterized by a operation point (Fi). The passage of 
the operation point from F 0 to F i is the transient state which 
gives place to fugitive primary effects. If the point of 
operation is unstable, the secondary effects appear as 
consequences of instabilities of the operation point. The 
following figure illustrates the transition from the system of a 
state F 0 to a state Fi. A defence plan can act on the causes by 
preventive measure or on the effects by curative measure. 



components to guarantee the service awaited and 
minimized the number of initiating events. 

S To guarantee a quasi-absolute permanence of certain 
vital functions even in the event of failure of the 
equipment which fills them. This is obtained by the 
material and functional redundancy. 

The reprogramming of the production constitutes the must 
securing task for the management of the electric system into 
preventive. The development of a new safe order moves the 
operation static point towards a state in which safety (N-l) 
and the operation margin are respected. This action is known 
as in closed loop since it controls the results. 
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Figure 2: Preventive action in closed loop 

2.2. Curative defence : this field gathers the whole of the 
automatic or manual actions which make it possible to detect 
the variations on certain sizes characteristic of the correct 
operation and to start the suitable corrective actions if 
necessary aiming at ensuring the protection of the materials 
and the safety of the system. The aim is above all to avoid 
the degeneration of incidents and/or failure in incident of 
great width. The counter measure adopted for these extreme 




Figure 1: Possibility of action defence plan 

2.1. Preventive defence: It consists into the actions carried 
out in order to avoid the starting of the collapse phenomena. 
The electric system must be insensitive to the consecutive 
losses of certain equipment following failures and/or risks 
considered to be probable and taken into account in the 
dimensioning of the system. 

The carried out actions aim mainly: 

S Preventive maintenance: to maintain the level of 
reliability, availability and performance of the 



Figure 3: Curative action in open loop 

There are 2 possibilities of curative action into curative, one 
based on the primary effects and the other on the side effects. 
However, the primary effects appeared in transient in the first 
moment and do not constitute a state of steady balance for 
the system. In the majority of the cases, these primary effects 
do not have consequences. Thus, they do not have direct 
impact on the safety of the electrical supply network and 
consequently they do not affect the customers. The action of 
such a plan based on the primary effects is known as in open 
loop. 

/Sh. 
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The disadvantages of this plan lie mainly in: 

• the difficulty in detecting the primary effects, 

• calculation must be very fast, 

• the risque of unnecessary operation: the primary 
effects may not have serious consequences (no 
appearance of side effects). 

The primary effects, in certain cases and in the presence of 
the worsening factors, can excite the electric system to give 
rise to side effects. The action of a defence plan based on the 
side effects is known as in closed loop. The control of the 
results is based on the calculation of the operation point as 
shown in the following figure. 




Figure 4: Curative action in closed loop 

3. Implementation of the defence plans 

The implementation of defence measures and the associated 
difficulties are exposed in this part while referring to the 
existing defence plan in the world integrating defence 
measures centralized and decentralized. 

We distinguish two groups of defensive measures: local 
protections and traditional which are generally recommended 
for the lines of interconnection and centralized defence 
measures. The following paragraphs introduce these two 
groups. 

3.1. Local measures of defence 

The implementation of these measures requires detailed 
specifications taking into account the specificity of each 
country and the already implemented defence plans. 

It is thus, to avoid the harmful consequences when a 
transmissible maximum is reached, it is recommended to 
equip each line with interconnection of the following 
protections: protection of maximum power, protection out of 
step and protection of under-frequency. The installation of a 
protection of low voltage is not recommanded, it could 
involve inopportune releases at the time of the transitory 
angular excursions. 

Protection of maximum power (maxP): The transmissible 
maximum power on an interconnection is the maximum 
capacity of exchange which can be exchanged between two 
zones without losing stability. Once this maximum is reached 
the system loses stability which results in a collapse of the 
voltage and/or a loss of synchronism. To avoid, in a 
preventive way, this type of break when a transmissible 
maximum is reached, a Watt-metric protection of maximum 
power (maxP) is necessary. 



The adjustment of protection maxP must ensure a sufficient 
margin compared to the maximum power transmissible. The 
threshold maxP is regulated below the computed value of 
the transmissible maximum. To take into account the 
influence of the topological changes and the variations of the 
voltage, a margin of 20% compared to the calculated 
maximum transfer is generally recommended. 

Out of step (max0): This protection must detect the 

conditions of loss of synchronism inter-zone or inter-area and 
separate them quickly. Its action is curative contrary to the 
protection maxP which has a preventive action (which acts 
before the transmissible maximum is reached). The 
protection usually used to detect the losses of synchronism is 
’’out of step protection " which is based on a measure of the 
impedance seen on the line at the some time detects, in the 
plan (Resistance: R, reactance: X), if a preset zone is crossed. 
Protection is activated when this crossing was carried out in 
a time higher than a minimum. The publication " Modeling of 
out of step conditions in Power systems " presents the 
principles and of the examples of application. 

Protection of under-frequency (minF): This protection 

ensures the separation of part of the load to restore in under 
network the balance Consommation/production. It limits the 
propagation of the frequency collapse from an area to 
another if a significant deficit of production is detected. The 
adjustment of this protection is generally located below a 
limited number of common thresholds solidarity (generally 
2) between the various inter-connected parts. 

3. 2. Centralized defence plans 

Coordinated measurements of defence are generally called 
special protection system (SPS). SPS is a protection system 
which was conceived to detect particular situations of the 
system which are known to cause an unusual stress of the 
system and which takes certain preset measurements to face 
the conditions observed in a controlled way. In certain cases, 
the SPS are conceived to detect operating conditions which 
were identified to cause an instability, an overload or a 
voltage collapse [ 27 

A synthesis of the SPS based on a collection of 
approximately information from 49 companies of 17 
different countries for a total of 1 1 1 SPS made it possible to 
raise the following strong and weak points of the these 
protection systems: 

S The rejection of production and the load shedding are 
the most current measures followed in the SPS; 

S The conditions of operating requiring the action of the 
SPS are not often met, but when it is the case, the SPS 
functions generally correctly; 

S Much SPS at any moment requires to be in service; 

S The alternatives which could be considered in the place 
of these SPS, are generally load shedding, the 
modification of the production and/or more restrictive 
operating condition. This last point means limitation of 
the capacities of transfer; 

S The cost of a disfunction of the SPS is generally largely 
lower than the cost generated by the situations than the 
SPS makes can avoid. This implies that the installation 
of SPS is economically interesting and profitable, taking 
account of the disfunctions; 
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V All the existing SPS are solutions dedicated for 
particular systems. No supplier offers a complete 
solution covering the problems of data acquisition to the 
execution of the defence actions while passing by the 
design. 

4. Defence Plan in the maghreb electrical system 

In this paragraph we present the existing defence plans of the 
maghreb inter-connected networks (Morocco, Algeria, 
Tunisia and Libya) as well as the advantages and 
insufficiencies shown by simulations and of experience 
feedback. In addition we studied the real recording of the 
active transit at the level of the Interconnection Morocco- 
Spain (IME) following the release of a group in Morocco and 
the curve of corresponding simulation after a series of 
adjustment of dynamic parameters. 

Releases of the Maghrebin interconnections during the 5 last 
years were presented. For the optimization of the loads to be 
lightened, works of simulation were started for a loss of 
production of 20%, 62%, 75 % and 80% of the peak power. 
The results show that a over shedding the two cases 62% and 
75% of the load is recorded. However, the use of a criterion 

df 

— based on the rate of frequency variation for the criterion 

dt 

of load shedding is recommended to optimize the load to be 
lightened according to the variations of the system 
parameters as well as the dynamics of the disturbance 
(quantity lost according to time). 

4.1 Adopted methodology of adjustments on installed 
protection of the maghreb defence plan 

In the electrical supply networks with strong interconnection, 
three categories of state are retained: 

V Normal operation: for which only actions of 

adjustment are taken (generally the frequency is in a 
beach of ±200mHz); 

V Operation in incident, which implies the use of the 
principle of solidarity. The adjustment of the frequency 
generally rests on the control of the production 
(generally the frequency is in a beach of ±lHz); 

V Operation in defence which relates to states out of the 
usual criteria and implements average extremes like the 
automatic load shedding and the isolation of the thermal 
units with their auxiliaries. 

In the case of network structures with simple interconnection, 
it is possible to establish a diagram of preventive shedding 
engaged simultaneously with the decoupling from the 
external network. Ideally, this initial shedding must be 
supplemented by a classical plan of shedding metric- 
frequency. 

In order to adopt the provisions to take within the framework 
of the defence plan as well as the adjustments to make, 
studies of operation of the electrical supply networks are 
carried out in order to release the constraints of the system 
following the great disturbances. The constraints are of two 
types: 

Static stresses: they are primarily the constraints of overload 
on the transport equipment or the voltage drops. 



Dynamic stresses: these constraints are related to the 

dynamic behaviour of the electrical supply network to the 
disturbances which results in three phenomena: 

Voltage Collapse 
Collapse of frequency 
Loss of synchronism 

The following figure shows the geographical diagram of the 
maghreb. 




Figure 5: Geographical diagram of the maghreb 



The operation studies of the maghreb network are carried 
out by the simulations with a function using the electric 
computation softwares requiring data bases of the various 
elements of the network (line, transformer, group... etc.). 
The dynamic behavior of the system depends on the 
response of the alternators and their systems of excitation 
and speed regulations, the parameters of the models used 

must be adjusted in order to obtain answers close to reality 

[ 28 ]. 

The approach adopted for this practical validation passed 
by the three 3rd following steps: 

1. To analyze the history of the incidents in order to 
draw real facts of the lesson being able to explain 
the variations observed between reality and theory; 

2. To make a study of sensitivity allowing to 
determine the most influential sizes following the 
loss of a group in an inter-connected system 

3. To check on real incidents (of which the data are 
available) the validity of the conclusions drawn 
before. At the time of this third phase, the most 
influential parameters are adjusted until obtaining 
results nearest possible to recordings. From the 
various simulated incidents, it is then a question of 
drawing the most adequate models to the maghreb 
electric system. 




Figure 6 : adjustement between reality and simulation 
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As an indication, the following figure represents a real 
recording of the active transit on the level of the 
interconnection Morocco-Spain (IME) following the 
release of a group in Morocco and the curve obtained by 
simulation after a series of adjustments of the dynamic 
parameters. 

4.2 Load shedding by minimum frequency (minF 
protection) 

Because of the electrical interconnection between maghreb 
and Europe (UCTE) the load shedding by relay with 
minimum frequency in the maghreb countries will be 
requested only in the event of significant deficit of 
production and the separation of the concerned network or 
the networks from UCTE system. 

In order to optimize the operation of the maghreb 
interconnections and to especially ensure a mutual help 
between the maghreb countries in disturbed mode, it was 
retained to adopt thresholds of load shedding to solidarity. 
This makes it possible to lighten the load on the whole of the 
maghreb countries inter-connected in order to adjust the 
consumption production balance and to preserve the maghreb 
interconnections and of the thresholds suitable for each 
network to safeguard the system in the event of deficit of 
production after the opening of these interconnections. The 
solidarity thresholds of load shedding and their lightened 
powers are as follows: 

S Thl: Fl=49.3 Hz: load shedding of 8% of the peak load . 
S Th2: F2= 49 Hz: load shedding of 8% of the peak load. 

S Th3: F3= 48.7 Hz: Opening of interconnections. 

Below the F3 (frequency of network separation), some 
thresholds are integrated in order to safeguard the overdrawn 
network (objet of the initial disturbance). These thresholds 
are studied according to the size of the means of production 
installed and the time of launching equivalent to each 
network. 

4.3 Load shedding controlled to the active transit on the 
interconnections and wattmetric defence plan (maxP) 

All losses of production on the maghreb system is 
completely (96 to 98%) filled by the European network. The 
help of the European network comes through the cable 
Morocco-Spain, it crosses the maghreb networks to the 
overdrawn country. The request of the maghreb 
interconnections is all the more strong when the lost power 
is large. Defence Plans (by wattmetric protection and loss of 
synchronism) are installed on the lines of interconnection and 
have the function of opening the lines of interconnection and 
limiting the incident and its propagation from one network to 
the other. After the wattmetric defence plan action, the 
network origin of the incident is overdrawn and the 
automatic plan of load shedding (by minimum frequency) 
will be put in action if the fixed thresholds are exceeded. 

In complement of the wattmetric defence plan, a plan of load 
shedding by control to the active transit on the lines of 
interconnection is installed with an aim to decrease the 
transit on the interconnection lines below the thresholds of 
the wattmetric defence plan and ensure the help and the 
connexity of the networks. This plan of shedding is fixed 
according to the thresholds wattmetric strown on the 



interconnection and the sizes of the groups installed (in order 
to ensure the safety N-l groups). 

4.4 Protection against the loss of synchronism (max0) 

Protections against the loss of synchronism installed at the 
level of the maghreb interconnections are of type DRS. The 
curve presents a recording of the voltage at the frontier station 
of the interconnection Morocco- Algeria following the release 
of two groups in Morocco at 26 january 1995 associated to 
releases of lines 225KV. 




Figure 7: voltage oscillation againt out of step 



The following figure represents a recording of the 225kV 
voltage at the of interconnection station Tunisia-Algeria 
following a prolonged short circuit on a line serving the 
frontier network on April 09 / 2001 . 




Figure 8: Volatge at the 225kV interconnection bus 



Current adjustments: 



Table 1 : 





SPAIN 


MOROCCO 


ALGERIA 


TUNISIA 


MOROCCO 


65%Un, 

2beat 


- 


70%Un, 

lbeat 


- 


ALGERIE 


- 


80%Un, 

lbeat 


- 


80%Un, 

lbeat 


TUNISIA 


- 


- 


80%Un, 

lbeat 


- 



4. 5 Load shedding by minimum valtage 

The load shedding by minimum voltage uses relays with 
minimum of tension and allows: 

S to avoid the electrical supply to customers under 
abnormal conditions of voltage. 

S to help the principal devices of the defence plan in 
failure. 

Moroccan network: 

- 1 5 stations are equipped by relays with minimum voltage to 
especially lighten local loads estimated at 20% of the peak 
power in the South network. These devices were set up in 
order : 
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S To especially avoid the collapse of the voltage at the 
level of the south network in the event of loss groups in 
western north 

S To ensure the help of the system of load shedding 
controlled Wattmetric protection on the IME. 

Algerian network: 

- 8 stations hight voltage and 12 stations middel volatge are 
equipped by relays with minimum voltage 

Tunisian network: 

- 2 stations are equipped with a relay wich regulated 
minimum voltage with 85 %Un and which lightens 40% of 
the load in each station. 

Libyan network: 

- Under study in 10 stations especially on the Libyan south 

4. 6 Release of the Maghreb interconnections: 

The table 2 gives the number of releases of the Maghreb 
interconnections by type of protection during the period 
2002-2006. 



of the load and the systems of voltage regulation and speed 
prove to be necessary for the study group of the Commission 
of the Maghreb Interconnections. 

5 Improvement activities and future recommendations 
5.1. Optimization of the load shedding by mixed criteria 
The power lightened in an electric system by a shedding by 
conventional minimum frequency (actual plan of shedding 
on the networks maghreb) is no optimal. It can lead to over/ 
under shedding of load sometimes nonacceptable. The use of 

df 

a criterion — based on the rate of frequency variation is 

dt 

recommended to optimize the load to be lightened according 
to the parameters variations of the system as well as the 
dynamics of the disturbance (quantity lost according to time). 

The following figure shows the results of simulation in the 
event of production loss of 20% of the peak load for various 
loading cases 80%, 75% and 62% of the peak load in the 
Tunisian network, in the event of use of a traditional plan of 
shedding frequency - metric: Over shedding noted in 
simulations were the cases 62% and 75%. 



Table 2: 
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4 


0 
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0 
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Figure 9 : Trip of tie lines interconnection in the electrical 
maghreb system 




Figure 10 : Mixed criteria of the load shedding 



The division of the shedding plan in two stages is necessary. 
These stages are: 

S A first part which uses the conventional shedding 
functioning in the event of loss of low power and 
moderate power; 

S Another which uses a mixed criterion of frequency and 
its derivative functioning in the event of loss of big size. 

The table 3 recapitulates the automatic load shedding plan 
suggested (result of the study) and compared to the 
recommendations of the UNIPEDE (Union of Production 
Energy and Distribution of Electricity). 



4.7 Coordination and exchange of information 

The design of the defence plans and the determination of the 
adequate adjustments are significant but not sufficient. 
Periodic controls of operation for devices of safeguard of the 
network as well as the sensitizing of the personnel in charge 
of these controls, are essential for the reliability of the 
electric systems. 

The defence plans of the grid systems require studies of 
operation and a control of the modeling of the inter- 
connected electric systems. The continuous effort for the 
development of experience sharing in the field of modeling 



5.2. Defence Measures centralized in the Maghreb 
interconnections 

The distribution of the adjustments on the various lines of the 
interconnections is a difficult spot since the distribution of 
the transits vary according to the production plan and the 
started production. The studies common at the maghreb scale 
(within the framework of the COMELEC) recommend a 
centralized measure of the sum of the transits on the various 
tie lines. 
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Table 3 : 



Thresholds 


Load shedding 
in % 


Recommendations 

UNIPEDE 


N° 


F(Hz) 


T(sec) 




1 st threshold: 
49 —48.8 


1 


49.30 


0.2 


5.4 


49.25 


0.2 


0.7 


Total of the threshold 


6.05 


2 


49 


0.2 


5.9 


Total of the threshold 


5.9 


15% 


3 


48.7 


0.2 


- 


2nd threshold: 
48.7 —48.4 


48.7 


10 


- 


4 


48.5 


0.2 


12.5 


Total of the threshold 


12.5 


15% 


5 


48.25 


0.2 


10.7 


3rd threshold: 
48.3 —48 


Total of the threshold 


10.7 


6 


48 


0.2 


9.3 


Total of the threshold 


9.3 


15% 


7 


47.75 


0.2 


11.9 




47.7 


0.2 


- 




Total of the threshold 


11.9 




Total load shedding 
(in %) 


56.25% 


40 —50% 



5.3. Future strategy of centralized defence 

If local measures of defence can be qualified as corrective 
and functioning in closed loop, centralized or cordinated 
measures can be qualified as preventive and generally 
function in open loop. Indeed, they aim at keeping the 
integrity of the inter-connected system acting quickly in the 
initially disturbed networks to ensure its stability. They are 
qualifieds as centralized because they imply the 
centralization of information towards one or more decision- 
making centres and thus of telecommunications. They are 
also characterized by the pre calculation of measures of 
defence according to the disturbance met. They will be 
necessarily doubled by local defence measures which fill in 
this diagram the role of back up protection. The presence of 
centralized measures of defence requires the temporization of 
the local protections in order to allow time to the system to 
react to centralized measures and give to local measures a 
role of safeguard. The diagram of safety of the system thus 
becomes organized in 3 levels. 

6. Conclusion 

In conclusion, we think that it is not possible to be denied 
against any type risk. Indeed, we cannot conceive all the 
combinations of breakdowns or incidents likely to intervene 
on such a large number of components on the one hand, and 
it is not economically justified to defend oneself (if we want 
to preserve a normal operation of the system) against risks of 
which the probability of occurrence becomes extremely weak 
and independent. Thus, for combinations of particularly 
severe risks but far less probable, we accepts degradations of 
the operation of the system leading to significant effects on 
the customers. The priority is then to preserve the control of 
the evolution of the incidents in order to limit their final 
impact. In more serious cases, one possibly agrees to 
sacrifice a reduced part of the system. 



The design of the Maghreb defence plans is based on the 
same techniques of safeguard, nevertheless the application of 
these techniques depends on the structure of the electric 
system: topology of the network, site of production. . . etc. 

The adjustments of protections adopted on the lines of the 
Maghreb interconnections, although they did not make it 
possible to ensure the mutual help and solidarity in term of 
frequency, showed a certain effectiveness to avoid the 
propagation of the incidents to the close countries 

With the development of the maghreb electric systems by the 
installation of the power stages of the significant groups, the 
reinforcement of the internal networks and the possibilities of 
extension of the synchronism zone, the defence plans must 
be revised and coordinated permanently and this in order to 
support trade and the mutual helps while preserving the 
safety of the electric systems. 



Table 4 : 



Level 


1 


2 


3 


Nature 


Preventive 

liability 


Preventive 

credit 


Corrective 
measure credit 


Objective 


Maintenance 
of the system 
inter- 
connected 


Maintenance of 
the system inter- 
connected 


The stability of 
the system inter- 
connected is lost. 

Opening of the 
sections critical to 
stop the losses of 
synchronism. 


Activation 


Disturbance 
less severe 
than the 
normative 
disturbance 
(called 

dimensioning 

incident). 


Disturbance 
more severe 
than normative 
disturbance. 
Activation can 
be anticipated 
for disturbances 
less severe that 
the disturbance 
normative. 


Disturbance more 
severe than 
normative 
disturbance and 
absence or faulty 
operation of 
measurements of 
centralized 
defences. No 
respect of the 
margins of safety 
on the critical 
sections, 


Means 


Definition of 
the incidents 
dimensioning 
the Margin on 
the capacities 
of transfer 


Tele shedding 
of load to 
remain below 
side levels of 
exchanges max 
on the sections. 
Fast decrease 
(transitory or 
permanent) of 
the produced 
mechanical 
power 


Opening of the 
sections on local 
criteria shedding 
of load Parcelling 
out of the network 
in zones Preset 
isolation of group 
of production 
from their 
auxiliaries 
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Abstract 

This paper is concerned with the control problem 
synthesis for nonlinear systems using linear technique. 
The proposed algorithm consists on an approximation 
of the non-linear systems by a set of uncertain linear 
models and presents a systematic approach for non- 
linear control design by using the gain scheduling 
technique to insure the transition of a non-linear 
dynamical process from an actual operating condition 
to a desired one. The non-linear system is represented 
in neighbourhood of equilibrium points by a family of 
uncertain linear systems. A robust control law is built 
so as to insure asymptotic stability to a given 
equilibrium within a maximal ellipsoidal region. Given 
a pre-specified equilibrium curve connecting the initial 
and final points. Output feedback for the local robust 
control synthesis are considered, together with local 
performance criteria. A simple numerical experiment is 
provided to illustrate both the effectiveness of the 
synthesis and the performance achieved. 

Keywords. Robust control, Gain Scheduling, Transition 
Control, Norm bounded uncertainties, Linear Matrix 
Inequalities (LMI). 

1. Introduction 

The majority of control laws implemented in industry 
are still designed and analysed using predominantly 
linear techniques applied to linearized models. Given 
the continuous increase in the demands on their 
performance and reliability, control law designers are 
highly motivated to explore the applicability of more 
powerful design and analysis tools allowing to take into 
account both nonlinearities and parametric variations. 
In this context, new robustness analysis and robust 
control synthesis methods are first developed for 
systems presenting an adaptive controller or a varying 
controller structure. The common approach to solve a 
complex problem is to use the strategy of divide and 
conquer. This strategy leads to multiple models or 
multiple controllers approaches to deal with nonlinear 



systems. Some of the approaches used frequently are: 
multiple models adaptive control, supervisory control, 
the quasi-LPV approach [19] and gain scheduling [18]. 

Gain scheduling dynamic controllers for a nonlinear 
plants are studied in [3]. One can find in [14] 
approaches to design controllers for linear parameter- 
varying systems. Takagi-Sugeno fuzzy gain-schedulers 
are proposed in [1] and in [12]. A generalized gain 
scheduling approach is developed in [6]. 

Gain scheduled control technique is based on the 
existence of several, generally linear, models which are 
assumed to be valid for the representation of a system 
to be controlled in different subsets of the state space. 
Obviously, for a complete characterization, the union of 
the subsets covers the state space domain of interest. 
For a nonlinear system, the models can be obtained by 
lineariation at several different feasible equilibrium 
points, referred to by trim points in the sequel. 
Different controllers are designed for each linear 
model. 

In this work, we develop an alternative gain scheduling 
procedure that exhibits some of the performance 
properties of standard gain scheduling, but provides 
stability guarantees between the scheduling points. In 
particular, we consider the use of the gain scheduling 
approach to the problem of stable guaranteed transition 
from an actual operating point to a desired one by using 
a pre-specified path in the equilibrium manifold of the 
non-linear system connecting the operating points. An 
uncertain linear system is used to describe the non- 
linear system in a ellipsoidal region around an 
equilibrium point. A Lyapunov function for the 
uncertain linear system is used to estimate ellipsoidal 
stability regions of the non-linear system. Repeating the 
procedure on the pre-specified path in the equilibrium 
manifold connecting the two operating conditions, we 
can cover the path by a set of ellipsoids. The algorithm 
developed assure closed loop stability of the non-linear 
system and local performance . 
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The paper is organized as follows: Section 2 gives a 
description of the HL-CPWL approach. Section 3 
concerns the local robust control laws synthesis : robust 

H 2 filter. We end the paper by an illustrative example 
and a conclusion. 



2. A High level Canonical Piecewise 
modelling: HL-CPWL approach 

Consider the class of discrete nonlinear systems 
described by the following state space equations and 
the following remarks: 

x = g (x (t),u (t)) + B (t) + B 2 u (t) 

y -C x {t) 

< 

X = [xj X 2 ..JC Jei?” U - [ft, U 2 --U m ](E R m 

y tR q w gR p 
with 

{g(x(t),u(t))=A n x(t) + B n h(x(t ),u (t )) ( 2 ) 

where 

h(x( t ), u (t) ) -> h(^(t))=[h, (m h 2 (m- h n m)f (3) 

g (x ) is a nonlinear function. 

x e R n , u e , y e and w e are respectively 

the state, the control, the measured outputs and the 
perturbation. 

We are interested in continuous time SISO or MIMO 
state space representation of systems. 

Systems are state nonlinear and/or control nonlinear. It 
can exist one or more non-linearities depending on one 
or two state variables and/or control variables. 

The linear part accepts a non null dynamical matrix 



The main steps of the procedure are described below 
and the next subsections are devoted to large 
explanations of each step of the HL-CPWL approach. 
The algorithm is composed of four steps 

2.1 Approximation of nonlinearities function 
In this section we introduce the basic aspects which are 
common to all types of PWL representations. We begin 
by briefly recall PWL functions, intersection and 
Simplex[7]. 



Definition 1 : PWL functions 

Let the domain is partitioned into a set of convex 
polyhedrons called regions E (0 so that 

S =U £ ' (0 =E m vjE {1) vjE (N) (4) 

i= 1 

by a set R = { R t a S J = 1...N ’ } of a finite number 



of -dimensional hyperplanes, also called boundaries: 



R f = { x ^IR n : 7r i (x) = aj x - /?. = 0 } (5) 

where a t eIR n , /?. elR Vz= 1.2V. Then, a PWL 
function is defined by the local (linear) functions : 

f^ l \x)=J^ l \ + cjh ) (6) 

withx tIR mxn and^7 (0 eIR m . Moreover 

J (,) jc + fl7 (i) =J (k) x + fi / k) ViEf (i) nE (k) 

if 

Definition 2: Intersection 

In a partition of the domain D(^R n ,ann-\ 
dimensional hyper-plane is said to be a first-order 
intersection, denoted by S (1) . A linear manifold of 
dimension (n-k) in D is a kth -order intersection 
S {k } if it is the intersection of two or more linear 
manifolds of type S {k ~ 1} , ie., 

S (k) :=f| Sf-v (7) 

2=1 

Definition 3 : Simplex 

Let x 0 ,x 1 ,...x n be n + 1 vectors in R n . A simplex 
S {l) = A(x 0 ,x l ,...,x n ) is defined as 

A(x 0 ,x 1 ,...,x„) = ^x eIR n :x =^/ 2 ,x,.| (8) 

n 

where <1 and =1, for i =0,1. ...ft 

2=0 

The interest of these HL-CPWL functions comes from 
their systematic writing for any functions. Then, it will 
well suit for functions, which represent non-linearities 
in systems. 

We consider rectangular compact domains in R n of 
the form 

S ={ x e/ 7T : 0 < x . < m i cr i ,z =1,2,..., ft } (9) 



Where <J i and m i are a grid step and the number of 
cells associated with x l axis. 

Each dimensional component x* of a compact domain 
S can be subdivided into m i subintervals of amplitude 



a i = (x + x)/ m i . As a result, S turns out to be 
partitioned into ]^[ i m i hyper rectangles and to 



contain N= Y\ m , +1 vertices. 

For as example, we introduce a non-linear function 
defined in R 2 . The simplicial partition of the hyper 
rectangle [0, o x ] x [0, cr 2 ] in the simplices 



A(0, (cr, , 0), (cr, , a 2 )), A(0, (0, a 2 ), (cr, , <r 2 )). 
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Figure. 1 A Simplicial partition of a domain in R 2 




It is important to note that the HL-CPWL function 
possesses the minimum number of parameters required 

to represent the PWL function^ . This follows from 
the fact that all the information that is necessary to 

characterize the function & is given by the function 
values on the vertices. 

With a fixed number cr. of hyperplanes, one can 

design the approximation g HLC pwL °f g , for an 
arbitrary precision that 

gHL -cpwl (x ) = C A(x ) (10) 

The matrices A(x ) and C T are constructed according 
to the HL-CPWL. [7, 15] 

Now, we can apply to the considered class of nonlinear 
systems (l)-(3) which is rewritten by 



Determination intermediate equilibriums points 

The intermediate equilibriums points will be 
determined as the middle of the validity domain of each 
hyperplans bounded by twoconsecutive breaks points 
BP j9 BP j+ 1- 

BP: + BP: +1 

X eq j= 1 2 7 (14) 

Dertermination the uncertainty domain 

The norm bounded uncertainty are defined as follows 

A A=D i F i E i (15) 



— S HL-CPWL (x) + B [W (t) + BjU(t ) 
1 y=Cx{t ) 



with g HL 

-CPWL (x)= 



1 HLCPWL 



kHLCPWL 



h 



\ nHLCPWL ) J 



(x) 

(x) 



(X) 



( 11 ) 



( 12 ) 



and h j HLCPWL (x ) A n x (t ) + b t (13) 

2.2 Design of uncertain LTI systems 
Let's consider the piecewise linear model (11)-(13) . In 
this section is given the way we have chosen the breaks 
and the intermediate equilibriums points from the 
piecewise-linear approximation of non-linearity’s. 

Determination breaks points 



In the subsection is given the way we have chosen the 
break ( BP j , BP- +1 . . . ) points from the HLCPWL 
approximation of the non-linearity. Figure 2 depicts a 
nonlinear function g(x ) , g hl-cpwl ( x ) an d break 



points BP { . 

These points are considered for each model uncertainty 
: limits validities domain. 



Di e TC* r and Ei e ${ l * n are the matrix constants. The 
uncertain matrix belongs to the following: 

F t ={F t eST*' :F?F t </} (16) 

The calculus of each ellipsoid domain of uncertainty of 
each non linearity is 



if g z (x) is one argument function x f =[x lz ]then 



D i =1 and E t = a t with a, = X eq i+ , ~X eqJ 
if g t (x ) is two argument function x = [x X i x 2i ] then 
D i =1 , and E t = [a i ,P i ] with Ot i and P t are 
represented respectively the two axis of ellipse ^ 



a . = x 



eq ,1 ,i +1 



" eq ,1 ,i 



and 



Pi ^ eq ,2,k +1 X eq ,2,k 

The last condition for the adjacent equilibriums points 
is mandatory for being able to apply the switching 
between two state space subsets. This condition is also 
important for the stabilization, which is described, in 
the following step. 

By this way, we can consider, for each state-space 
subset, one equilibrium point and associate one 
uncertain local linear model 



\x(t) = (A .+A4)x(0+ B 2 u ( t) +B x w (t ) +b t 
{ y (t)=Cx(t) 



V/=1..JV(17) 



Vx(t) e X f 

With Xj is the variation domain of the state variables. 
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The equation (17) represents a problem at the level of 
the term J3- . For this reason we must eliminate this 
term, in order to relocate the state vector x(t) by a 
constant xq . 

Let’s suppose now that, we have a new state 
variable z{t) where z(t) = x(t) - x Q , z(t) = x(t ) . 

Now, we can apply to the considered class of linear 
model (17) , which is rewritten by 

\z(t) = (A i +M i )z(t)+B 2 u(t)+B x w (t) 

{ y(t)=Cz(t ) (18) 

Vz(?) eZ ; - 

withZ z - the variation domain of the state variables 



2.3 Local robust synthesis 

The next section is devoted to large explanations of 
robust synthesis for a family of continues time bounded 
norm uncertainly linear systems. 

2.4 Global control law by scheduling 

The last step for this algorithm is to calculate, for each 
uncetain linear model, a constant local control law 
which is differents form on equilibruim point to 
anothers one. This family of constant local control laws 
permits to reach every equilibruim points from each 
adjacent equilibruim points taken as initial condition. It 
may take a long time to reach them, but we assume that 
when the current state is sufficient near the itermediate 
equilibruim point, the switching between two state 
space subsets may happen. By this way, the proposed 
approach ensures a robust control synthesis. The 
smooth swicthing is obtained with the sufficient 
conditions of inclusion. 



3. Robust H 2 Filter 

Let the uncertain continuous-time invariant linear 
system be described by 

x =(A f +D i FE i )x +B x w 

< y -Cx F'F<I (19) 

z =Hx 



where where, x ^R n , u e , y e and z e Tf are 

the state, the control, the measured output, and the 
output estimates vectors. Furthermore, w represents an 



exogenous signal and H zw (s ) denotes the transfer 
matrix function from w to z. 

The basic problem to be formulated is the stabilizability 
problem with respect to the strictly proper filter 



= A 0 x 0 +B 0 y 
z o = C 0 x o 



( 20 ) 



We can determine the matrices A Q ,B 0 and C Q if 
z Q ^z .For this reason we have identified a problem 
of type H 2 on the system extended 
x = (A +DFE)x +B x w 

< x 0 =A 0 x 0 +B 0 y (21) 

£ - FLx -C 0 x 0 

with £ the error signal and CO perturbation . 

In this section, a performance specification is added to 
the previous stabilizability problem, the problem being 

to find a filter of class (20) that keeps the H 2 norm of 

the transfer function H zw (s ) as small as possible. 

Provided matrix A g in (23) is Hurwitz, the square of 

this norm can be expressed in terms of the solution to a 
Lyapunov equation (observability Grammian) such that 
this minimization problem with respect to the triplet of 

variable (A 0 ,B 0 ,C 0 ) is given by 

mm(A 0 ,B 0 ,C 0 )ju 

VF : F'F < 7,|C g (pi -A g -D g FE g T l B g \^ <p (22) 



where 
A„ = 



A 


0 ' 


= 


V 


_B 0 C 


1 

o 


’ s 


0 



A, =[H -C Q ], 



° g = 



(23) 



and E = [E 0] 



The command with observers increased returns to the 
terms of this uncertainty as: 

min trace W 



W>B g PB g 

< p g + V, +C X +p g D A p g +E ' g p g <° 



(25) 



which, using Schur complement and transformation can 
be rewritten in the equivalent form 

wmtraceW 

(26) 



>0 



" W 


F 


B\X 


F 


Y 


I 


XB, 


I 


X 



AY +YA ' 


Z+M 


YH'-L' 


D 


YE' 


Z'+M' 


AX +XA+FC+C'F' 


H' 


XD 


E ' 


HY -L 


H 


-I 


0 


0 


D } 


DX 


0 


-I 


0 


EY 


E 


0 


0 


-I 



with 

L =C Q V ? , F = UB 0 , M -V A 0 U' and 
Z =A +YA X +YC F 



A numerical example to valid our method is provided in 
the next section. 
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4. Numerical Example 

Consider the Tree Tanks system shown in Fig.3, and 
described by the set of state space equation (27) . 






T \ h h 

n , 




Figure 3: Tree Tank system 



X 1 =“^13#1 ~ X 3\+J U l 



x 2 _ ^32\/| x 3 x 2 1 ^20V X 2 + ~ w 2 ^7) 



X 3 _ ^13\/ X 1 x 3 ^32\/ x 3 X 



with 



K 13 = 



^3 xS c xsgn (H l ( t)-H 3 (t))x^2xg 



K 32 = 



K 20 = 



^32 xS c X s g n (#3 (0 - H 2 (0) x V 2 *.? 



and 



^20 X S c X ^/2>< 



& 









x \ 


x= 


H 2 


= 


x 2 








_*3_ 



"a" 




u x 


_Qi_ 




u 2 



, w : 



where, x eR n and u e are respectively the state 
vector and the input vector. 

The system present two nonlinear functions two 
arguments and one nonlinear function one arguments , 
we consider the index r a number of nonlinear function 
in the system. ( r = 3 in this system ) . 

where D k = / 



and E k = 



with 






0 



0 ( a k + Pk ) 






Pk 

a k 

(«* +A) 



k 




CO 


a 


p 


1 


-0.774 


0.6754 


0.25 


0.62 


2 


0.5477 


0.368 


0.23 


0.53 


3 


-0.5476 


0.368 


0.29 


0.59 


4 


0.773 


06755 


0.32 


0.67 



For this example, we chose to limit the variation 
domain in levels of tanks at the interval ]0, 60] and the 

number of local models uncertain is 2 models (N = 2). 
The application of the HL-CPWL representation on the 
nonlinear functions 

h[ = Jx 2 i ,r = 1 ,/ = 1...2 , 



h [ = Jiy -x 3i \,r =2,i =1...2 and 



h i = y| X 3,/ - x 2,i\’ r=3 ’ i =1- 2 

Can determine the intermediate equilibrium points 
(table 1). 

Table 1.1 : Equilibrium points value 



i 


1 


2 


X eq ,i 


12,52 


37,52 



Fig.4(a -b) depicts the both non linearities h[ and its 
approximation. 

In this example we have 4 models in the following 
domain which corresponds to a good enough 
approximation of the non linearity. 

The analytic function of the nonlinear functions 
approximation is 

4 

/ hl-cpwl( x ) = Yj^k x i +t> k ) + (0) k Xj +S k ) i ={1,2 }and j ={l,2}‘ 
k= l 

The determination of different equilibrium points and 
the number of linear models can define the parameters 
of the uncertainly domain elliptical for each linear 
model. 

The uncertain linear multi-model constructed is given 
by 



(-K 13% 


0 


~ K u (°k " 




\ 


0 


-(K 32 G>k + E 20 7 1k) 


~ E 12 Vk 


+ D k F k E k 


Z 


[ ^13% 


Ky2 m k 


Ki3<n k -K 32 rjk y 







— o 

S' 



0 1 
s 

0 0 



Robust H 2 Filter 

It is interesting now to apply the proposed method to 
design the robust H 2 filter (26) for (27) formulated 
from (19) with a disturbance w and for the same 
nonlinearties. 

The control objectif is to drive the state variables 
\^H X H 2 H 3 J from an initial values [040 30] r to the 

target equilibruim points [29 19 30] 7 . 

Fig. 6a shows the switching between state subpaces and 
the both responses have reasonably good transient 
performances thanks to the uncertain domain choice 
and the Fig. 6b shows the gain scheduling control. 
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This technique ensures stable transition between the 
starting point and destination by four corrector. For the 
latter, there are discontinuities in commutation 
moments, the amplitude of these discontinuities 
decreases with the increase in the number of corrector. 
The algorithm permits to tune the complexity of the 
global control law and the gain scheduling technique 
with the cardinality of the set of local linear systems. 

5. Conclusion 

In this paper, we have proposed an algorithm that 
approximate non-linear systems by a high level 
canonical piecewise linear and the algorithm guarantees 
the closed loop stable transition of a nonlinear system 
from an actual operating point to a desired one. A norm 
bounded uncertain linear system is used to describe the 
nonlinear system in a ellipsoid region around an 
equilibrium point. The robust control laws scheduling is 
executed when the state trajectory reaches the attraction 
domain of the following equilibrium point. The 
scheduling policy is based on the belonging of the 
estimated state to the stability ellipsoid of the next 
equilibrium point. 

Although we have restricted ourselves here to a simple 
class of nonlinear systems, nonlinear systems with one 
non-linearity and two non-lineartity, such an approach 
can be extended to larger class of nonlinear systems. 
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Abstract 

Fuzzy modelling has been widely applied as a powerful 
methodology for the identification of nonlinear system 
from process measurements. This work proposes a 
contribution on the application of fuzzy logic to 
identification of a class of large scale system with 
multiple models. The researched models have black box 
structure and the rules are determined by numerical 
approaches. The proposed approach is developed by 
using fuzzy clustering and is applied on a real system 
which is nonlinear in nature (greenhouse). 

Keywords: Identification of nonlinear systems, Fuzzy 
clustering, greenhouse climate modelling. 

1. Introduction 

Recently, there is a huge demand for tools that can be 
used for identification of complex nonlinear process 
which are difficult to describe mathematically. In the 
literature technique for mathematical modelling for real 
process are classified in two categories: physical 
modelling and system identification [6]. The first is based 
in terms of the physical laws involved in the process and 
the second is based on the analysis of the process input 
output data and from empirical expertise by using fuzzy 
logic which effective for approximation of uncertain non 
linear dynamics systems. Generally, modelling process 
consists to obtain a parametric model with the same 
dynamic behaviour of the real process. However, when 
the process is complex, it is very difficult to define the 
different mathematical or physical laws which describe 
here behaviour [2], [1]. Numerous approaches have been 
proposed [4] which compute nonlinear dynamic fuzzy 
models from input/output measurement data [7] or neuro- 
fuzzy approaches [5]. 

The fundamental goal of this work is addressed of fuzzy 
modelling by using an algorithm that implement fuzzy 
clustering. 

The proposed algorithm is based on decomposition of the 
fuzzy relation into sub relations, through a process of 
unfolding the fuzzy rules. The obtained sub rules are then 
grouped into subgroups (clusters). The number of the 



clusters is defined on line. The clusters obtained with this 
fuzzy clustering algorithm are used to compose the new 
fuzzy submodels. The proposed approach is applied to 
model the greenhouse air temperature and humidity 
under process control. Each fuzzy submodel defined 
reflects the behaviour of one real physical sub model. 

The paper is organized as follows: section 2 explains the 
identification of fuzzy model, in section 3, we describe 
the greenhouse climate model, in section 4 we introduce 
a fuzzy clustering algorithm, simulation results are given 
to demonstrate the effectiveness of the proposed 
approach in section 5. 

2. Identification of fuzzy model 

Process modelling is a problem particularly significant 
for the implementation and the development of controlled 
systems. In this section it is assumed that fuzzy systems 
are multi input single output Y : U —> V where the 

input space is U = u l u 2 u N and the output space 

is V. 

Consider a system which is described in the discrete time 
by y(k + 1) = f(z(k)) (1) 

where the regression vector 

z(k) = \u x (k),...u r (k),q l (k\...q p (k\y x (k\...y m (k)\ (2) 

with r process input U, p measurable disturbance q, and 
m process inputs y. the unknown function f(.) is 
approximated by using Takagi-Sugeno type fuzzy model, 
represented by a collection of fuzzy if-then rules: 

R i : if zj is A n and ... if z nz is A hlz ^ 

then y, (k + 1) = a Uo + a m z x + --- + a Unz z nz 
with nz = r + m + p and A u inU are linguistic terms 
characterized by fuzzy membership functions (z) . 

1 is the index of rules (1=1,2, ,M) 



A/( z ) - ex p(- 



-1 

2 



(Zl+ ! n)2 ) + ... + exp(— (z ” z + 2" z)2 

a n 2 L 



( 4 ) 
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where ju J (.) is the degree of fulfilment of rules 1, the 

centers Cj and the standard derivations <J l are the 
parameters of the Gaussian membership function. 

The output of TS model will be calculated with the 
following equation. [ zeng 1994] 

M 

'Y J /Ui(z(k))y l (k + \) 

y(k + J ) = ~ Jf, ( 5 ) 

/= i 

M 

or y(k + \) = '^ ( p l (z,c l ,a l )y l (k + 1) 

1=1 



where ^(z,^/,^) is the validation function for the 
Gaussian membership functions with centers Cj and 
standard derivations cr/ defined as 



&l(z,Ci,l 7,) 



M Z W) 

M 

^/r,(2(/c)) 

i = 1 



( 6 ) 



3. Greenhouse climate model 

The greenhouse climate model describes the dynamic 
behaviour of the stated variables by use of differential 
equations for the air temperature, humidity, C02 
concentration etc. these differential equation are results 
of combination of the various physical processes 
involving heat and mass transfer tacking place in the 
greenhouse and from the greenhouse to the outside air 
Figure. 1. a:. 




Figure l.a: Scheme of greenhouse climate model 



The important number of variables, (references, 
perturbations and commands) and the complexity of 
phenomena (biologic, weather, evolution of the 
plants....), make that we have a system multivariable, 
nonlinear and non stationary. Moreover, the 
perturbations, as the wind velocity and the global 
radiation, can sometimes have a power more raised than 
the command such the heating. 




For these reasons, we thus preferred to apply a 
multimodal based on fuzzy logic and taken into account 
all the variables of which we dispose Figure, l.b: 
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Figure l.b: Input/output diagram of a greenhouse 



where , Ti and Hi are respectively temperature and 
relative humidity of the internal air, the perturbation 
variables are Te (external temperature), He (external 
humidity), Rg (solar radiation), Vr (wind velocity) and 
the input variable are Ch (heating), Br (moistening) and 
Ov (roofing) 

To simplify the complexity of the considered system, 
greenhouse, we have proposed to decompose the MIMO 
system into a MISO. By this decomposition, the MISO 
systems are represented by compact TS model and 
convert a nonlinear problem into linear one. 

The TS model is defined by the rules which are linear 
conclusions of the systems outputs, for example for the 
rule j: 

if o v (k),ch{k),...R d (k)is V. then T {k + \) 

( 7 ) 

= a jX o v {k) + a j2 ( k)ch(k ) H h a j9 R d ( k ) 

f o v (k),ch(k),...R d (k)is w, then H n (k + 1) 

= b n ° v (*) + b n ( k)ch(k ) + --- + h iq R d ( k ) 



with 

A A 

0 X = [a jX ,a j2 (k),...a j9 ] and d 2 =[b n ,b l2 (k),...,b l9 ] 

are respectively model parameters for the temperature 
and humidity for the rule j. V j and W l are membership 

functions centers which are defined later. In this work we 
are interesting on a nonlinear with time variant 
behaviour. There for, an on line adaptation is necessary 
to obtain a good model able to describe the process in a 
large operating points to be used in a schema of 
adaptation control [3]. 

The rule conclusion is a linear model (3) calculated 
through a last squares approximation method, 
considering the input data and the associated output. The 
estimated parameters are chosen such they minimize the 
criteria: 



J(0,(k)) = 

Y J h i {k)(y i -z T eyiy.-z 7 9 t ) 



where y t is the output 



and h t (k) = 



Ri(z(k)) 



M 



(< = 1,2 ). 



/= 1 



( 9 ) 
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// z is the membership function of the data Z(k) at the 
rule j. 

We put H t = diag (h { (k), , h R (k)) (10) 



A A 

The parameter vector 0 X and 0 2 are updated online, by 
consideration that only the rules consequence are adapted 
for each local model by a recursive of the weighted least 
squares algorithm with forgetting factor X : 



A A 



d i (k) = e i (k-\) + 8 i {k)(y i {k) 




A 


(11) 


-z T {k)9 i {k-l)) 




X _ p t (k-])z(k) 

O i T n 


(12) 


Z {k) P,(k-X)z(k) + — r 




Pi (k) 




Pl =j[I-S,z T (k)]p,(k- 1) 


(13) 



where h t (.) is the weighting of actual data with the rule 
activation and p f is a matrix of the adaptation gain 



4. Fuzzy clustering algorithm 

The objective of clustering methods is to perform a 
partition of the collection of elements into c data sets, 
where the number c of clusters is defined online and 
depend of the state of day. In the generality of the 
situations, these elements are vectors of points in a space 
with a norm and distance measured and its results is a set 
of cluster of discrete points. 

The objective of clustering is to partition a data set Z into 
c clusters. Through clustering, the fuzzy partition vector 
ju(z(k)) = [jUj(z(k))\ cxl is obtained, whose elements 

]U k represent the degree of membership of the 
observation z(k) in the cluster v y ( j - 1,2,..., c ) . 

The partition should find a vector of clusters center, the 
number of clusters and a partition vector ju(z(k )) . This 
procedure is composed into four steps: 



Step 1. Initialization 

For a set of point Z, we calculate the distance between 
z(i) and z(j). 

D 1 = ( z(i)-z(j) T (z(i)-z(j )) (14) 

we determine the data point which done the maximum 
distance. 

Dmax =, ma \X D y) 

(15) 

= (z(0 - z{j) T (z(i) - z(j)) 
where N is the number of data. 

we initialize the centers of the two first clusters v x and 
v 2 with V x = z(i) , V 2 = z(j) and C — 2 
Step 2. compare 
max(|L z - - z(k ) ||) to £ 

If max (||F Z - - z(&)||) > s c — c + i and V c — z(k) 



Otherwise compute the c, mean vectors. 

Step 3. compute the new partition vector jU and calculate 
6\ and 62 

Step 4. compute the output of system, go to step 2 for 
future measure 

In the simulation section, this strategy is applied to break 
the greenhouse fuzzy climate temperature and humidity. 

5. Simulation results 

In this section the results achieved with the proposed 
greenhouse climate strategies are presented. The 
identification of both models was done by using input- 
output data previously collected. 

In figure 2 and 3, we present the results obtained from a 
simulation compared during a test period. These curves 
show a good agreement between measured data and the 
models outputs. 

Figure 4 shows the evolution of parameters tetal and 
teta2 

6. Conclusion 

This paper proposes a study of the application of fuzzy 
clustering to the identification of Takagi-Sugeno (TS) 
fuzzy models. These local models correspond to the 
different rules generated automatically which have 
variable consequent linear parameters in contrast to the 
common approach using fixed consequent linear 
parameters. 

The performance of the proposed modelling technique 
was demonstrated and successfully applied to real 
problems such as modelling of physical and biological 
process. The fuzzy clustering has been tested to split the 
inside greenhouse air temperature and humidity. 
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Figure 2. Comparison of the process output and the fuzzy model output with parameter’s adaptation 
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Figure 3. Comparison of the process output and the fuzzy model output with parameter’s adaptation, (by considering consign 

T ic =20°c, and H ic =70% ) 
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Abstract 

This paper presents a practical modeling and 
solution approach for the cyclic inventory routing 
under constant customer demand rates. The 
solution approach presented in this paper assumes a 
given vehicle type. The heuristic solution approach 
is well capable of finding the appropriate cost 
trade-off under varying circumstances. Furthermore, 
when fleet sizing is not considered, the heuristic 
outperforms an existing heuristic in finding the two 
way trade-off between distribution and inventory 
costs. The objective of this inventory routing 
problem (IRP) is to determine a distribution plan 
that minimizes average distribution and inventory 
costs without causing any stock-out at the 
customers. Deterministic constant customer demand 
rates are assumed and therefore, a long-term 
cyclical approach is adopted, integrating fleet sizing, 
vehicle routing, and inventory management. Further, 
realistic side-constraints such as limited storage 
capacities, driving time restrictions and constant 
replenishment intervals are taken into account. 
Computational experiments show that our column 
generation based heuristic solution approach is well 
capable of finding the appropriate cost trade-off 
under varying circumstances. 

Keywords: IRP, VMI, Supply chain, distribution 
and inventory costs 

1. Introduction 

The problem of integrating inventory and distribution 
decisions has been approached in different ways 
depending on inventory policies used at the 
customers, service level restrictions, and the time 
horizon considered. Integrating decisions related to 
different stages allows trade-offs to be examined and 
a better overall performance to be achieved. Over the 
last decades, both academics and practitioners have 
become aware of the importance of collaboration 
between the different stages in a supply chain. 
Integrating decisions related to different stages 
allows trade-offs to be examined and a better overall 



performance to be achieved. One particular example 
of supply chain collaboration is the concept of vendor 
managed inventory (VMI), encountered in 
distribution planning. Under a VMI-strategy, 
customers make their inventory data available to the 
distributor, and leave the responsibility of deciding 
when and in what quantities they are replenished to 
the distributor. The distributor thus has more freedom 
in designing efficient vehicle routes. Benjaafar, Saif 
[1] determine an optimal demand allocation and 
optimal inventory levels at each location so that the 
sum of transportation, inventory, and backorder costs 
is minimized. Chen, Yee Ming and Lin, Chun-Ta [2] 
solve the multi-product multi-period inventory 
routing problem with stochastic demand. 
Computational results demonstrate the importance of 
this approach not only to risk-averse decision makers, 
but also to maximize the net present value at an 
acceptable service level. As a result, an optimal 
portfolio (R, s, S) system of product group can be 
generated to maximize the net present value under an 
acceptable service level in a given planning horizon. 
These seminal papers have incited a body of literature 
on inventory routing, most of them considering 
stochastic demands and short time horizons. Sun, 
Mao and Li, Wen-Nian [3] give An inventory routing 
problem for a two-echelon distribution system 
consisted of one supplier (or distribution center) and 
many retailers based on deterministic demand under 
the mode of traditional retailer management 
inventory was studied. Raa, Birger and Aghezzaf, 
El-Houssaine [4] discuss the challenging optimization 
problem that arises in this context, known as the 
inventory routing problem . The objective of this IRP 
problem is to determine a distribution plan that 
minimizes average distribution and inventory costs 
without causing any stock-out at the customers. 
Deterministic constant customer demand rates are 
assumed and therefore, a long-term cyclical approach 
is adopted, integrating fleet sizing, vehicle routing, 
and inventory management. Savelsbergh, Martin and 
Song, Jin-Hwa [5] develop an integer programming 
based optimization algorithm capable of solving 
small to medium size instances. This optimization 
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algorithm is embedded in local search procedure to 
improve solutions produced by a randomized greedy 
heuristic. We demonstrate the effectiveness of this 
approach in an extensive computational study. 

The cyclic inventory routing problem has to be 
analyzed from a long-term perspective, in which the 
vehicle fleet is variable. Therefore, decisions on the 
vehicle fleet size have to be incorporated and a fixed 
vehicle cost has to be taken into account. However, 
most of the above mentioned papers use single routes 
per vehicle and assume that enough vehicles are 
available. In other words, vehicle fleet sizing 
decisions are not considered and the assignment of 
routes to vehicles is not taken into account. 

Raa, Birger and Aghezzaf, El-Houssaine [6] proposes 
a practical solution approach for the challenging 
optimization problem of minimizing overall costs in 
an integrated distribution and inventory control 
system. Constant customer demand rates are assumed 
and therefore a long-term, cyclic planning approach 
is adopted. 

Zachariadis, Emmanouil E and Tarantilis, Christos D 
[7]. studies an inventory routing model which 
integrates two important components of the supply 
chain: transportation logistics and inventory control. 
The aim of the problem is to determine the timing 
and size of the replenishment services together with 
the vehicle routes, so that the total transportation and 
inventory holding cost of the system is minimized. In 
methodological terms, they propose a solution 
approach applying two innovative local search 
operators for jointly dealing with the inventory and 
routing aspects of the examined problem, and Tabu 
Search for further reducing the transportation costs. 
The proposed algorithmic framework was tested on a 
set of new benchmark instances of various scales. It 
produced satisfactory results both in terms of 
effectiveness and robustness. 

The capacity constraints result in a restriction on the 
time between consecutive replenishments of a 
customer, i.e. in a maximal cycle time. On the other 
hand, there is also a minimal cycle time, a minimal 
amount of time between consecutive deliveries 
because the vehicle needs some time to drive to other 
customers and to the depot before it can come back. 
However, this minimal cycle time is not found in any 
of the above mentioned papers. 

Shen, Sheng-Yuan; Honda, Masakazu [8] focus on 
developing an integrated replenishment and routing 
plan that takes into account lateral transfers of both 
vehicles and inventory for a three-echelon supply 
chain system including a single plant, multiple 
distribution centers and multiple retailers. A mixed 
integer programming model to the overall system is 
formulated first, and then an optimization-based 
heuristic consisting of three major components is 
proposed. Aghezzaf, E.-H [9] discusses the case 
where customer demand rates and travel times are 
stochastic but stationary, and proposes a version of 
the inventory routing optimization model that 



generates optimal robust distribution plans. The 
approach proposed to obtain and deploy these robust 
plans combines optimization and Monte Carlo 
simulation. Optimization is used to determine the 
robust distribution plan and simulation is used to 
fine-tune the plan's critical parameters such as 
replenishment cycle times and safety stock levels. 

The problem at hand is thus the cyclic inventory 
routing problem (CIRP), in which a vehicle fleet and 
a cyclic distribution scheme for the selected vehicles 
have to be designed such that the total cost rate is 
minimized. This paper extends earlier work [8,9] by 
presenting a practical modeling and solution 
approach for cyclic inventory routing problems that 
takes into account driving time restrictions. The 
remainder of this paper is organized as follows. 

Section 2 presents the modeling approach. Section 3 
describes the solution approach. Section 4 gives 
Comparison to another heuristic. Section 5 provides 
conclusions and directions for future research. 

2.Modeling Approach 

Under the VMI-strategy, the distributor is faced with 
the problem of integrating its vehicle routing 
decisions with the inventory management of the 
customers. This integrated problem is known in 
literature as the inventory routing problem (IRP). 

In the long-term, the vehicle fleet is variable and fleet 
sizing has to be integrated with routing. Therefore, 
the routing concept is generalized such that a single 
vehicle can make multiple tours with different 
frequencies. Solving the cyclic inventory routing 
problem with constant demand rates is highly 
complex. It can be considered as a combination of the 
following inter-related and sub problems: 

1) Partitioning and assigning customers to a 
number of vehicles. 

2) Further partitioning of the subset of customers, 
assigned to a vehicle, into sub-subsets that will 
each be replenished in a separate tour. 

3) Constructing for every customer sub subset the 
shortest depot-to-depot tour that visits all 
customers in the sub subset. 

4) Determining for each vehicle an appropriate 
cycle time and frequencies for its tours. 

5) Constructing for each distribution pattern a 
delivery schedule to check feasibility and deter- 
mine the cost rate. 

Consider a vehicle replenishing a subset of customers 
S. The vehicle makes a distribution pattern with a 
cycle time of T days, consisting of n different tours, 
each visiting one of the disjoint customer subsets 
S i (i £ 1,..., 72 ) • For each tour /, a (relative) frequency 

f. is determined, meaning that the tour is madeyj 

times during one cycle. Since the number of days 
between consecutive deliveries to a customer has to 
be constant, the number of days between consecutive 
executions of a tour also has to be constant. 
Furthermore, the driving time is restricted to H hours 
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per day. As a result, the number of days in one cycle, 
i.e. the cycle time T of the distribution pattern, has to 
be an integer multiple of all tour frequencies f . The 

least common multiple of all f. thus defines a base 
cycle time of D days. 

One execution of tour i takes 7" / . hours for 

making the TSP-tour through the set s+ , plus 
Yjjes+tj hours for delays at all locations in s+ , 
where 7- stands for the unloading time at customer 



j e s. , and Note that this tour duration cannot exceed 

H hours. A distribution pattern cannot be restarted 
before the vehicle has made f. times each of its 

tours i. Knowing that only H hours of driving are 
possible per day, and that the number of days in a 
cycle has to be an integer multiple of D days, this 

leads to the following lower bound T min for the 
distribution pattern cycle time T : 

1 



T =D 

mm 



HD 4 



I ^(W) + S*;) 






a) 



The maximal quantity that the vehicle can deliver to 
the customers in the tour is thus f. times the vehicle 

capacity C. To avoid stock-outs, the vehicle has to 
return to its customers before they have consumed 
this hill truckload. Therefore, T/ f. , the time between 

consecutive deliveries, cannot be larger than the 
vehicle capacity divided by the cumulative demand 
rate of all customers in the tour: T- < — , with 

Ji > d , 

La jeSj J 

d t expressed in units per day. Similarly, the limited 
storage capacity of the customer /, denoted by • 
Taking into account the base cycle time of D days, 
the upper bound T max on the distribution pattern 
cycle time T is as follows: 



T =D 



^-min( 

D iGl - n / 



Cf t • ftfj 
1 — ,mm — - 

d. J gS d. 

f jeS J J 



( 2 ) 



Obviously, a distribution pattern is only feasible if the 
lower bound is not greater than the upper bound on 

the cycle time: T mi <T miix . 

The cost rate R of a distribution pattern consists of 
the following components: 



1) 

2 ) 



The fixed vehicle cost of a $ per day. 

The variable transportation cost of 



l7^ ( ,-) $perday - 



TSP(Sl) 

1=1 

3) The fixed vehicle loading and dispatching cost 

n 

of ^ f i fi A $ P er cycle. 

i=i 

4) The fixed unloading cost at the customers of 



fPj $ per cycle. 

z=i 

The stock holding cost at the customers. The quantity 



delivered to customer j e S. ? denoted by q . , covers 



demand for -Jr days, so n . =d^- 

1 J J Ci 



The total cost rate of the distribution pattern 
replenishing the set of customers S is given in the 
following formula: 






(3) 






ypj 



JjsS ‘ 2/,. 

This cost rate varies with the cycle time and with the 
frequencies f. . Thus, for every combination of 



frequencies, there is a theoretical optimal cycle time, 
for which the cost rate is minimal. This occurs when 
holding costs at the customers are balanced with the 
sum of transportation costs and fixed tour costs 
(components (2)-(4)). This is a generalization of the 
well-known Economic Order Quantity model, and 
therefore, this optimal cycle time is called the EOQ 
cycle time and denoted T EO q: 



T = 

EOQ 



Z /(V, + EkA) 



zz 



rA i 

‘ itS ‘ 2f, 



(4) 



Once frequencies f. and a cycle time T of a certain 



number of days are selected, each execution of every 
tour in the distribution pattern is assigned to a 
particular day in the cycle such that (1) the number of 
days between consecutive executions of the same tour 
is constant, (2) the vehicle drives at most H hours per 
day. Because of these two constraints, it is unlikely 
that a vehicle drives around the whole time and thus 
idle time needs to be introduced in the vehicle’s 
cyclic schedule. 



3. Solution Approach 

As the modeling approach presented in Section 2, this 
section discusses a solution approach that tries to find 
good-quality solutions that achieve good trade-offs 
between the different cost components. The solution 
approach presented here extends previous work [8,9] by 
considering the more general routing concept of 
distribution patterns and taking into account driving 
time restrictions. 

Here, an improvement heuristic is proposed that can 
be applied to an existing solution at any time to 
improve the solution quality. The improvement 
heuristic consists of removing and re-inserting a 
customer. Re-inserting a customer into the solution is 
similar to inserting a customer in the insertion 
heuristic. The customer is inserted into the existing 
distribution patterns of the solution, both in a separate 
tour and in the existing tours, and the cheapest 
alternative is kept. If a customer ends up in the same 
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position as before, the solution is restored and no 
improvement is found. If the customer ends up in a 
different position, the solution has improved. The 
order in which the customers are considered for 
removal and re-insertion is also determined by the 
customer priority. 

We need to find a collection of distribution patterns, a 
partition of the set of customers S into disjoint 
subsets, each covered by a separate distribution 
pattern, needs to be found, such that the sum of cost 
rates of the distribution patterns is minimal. If we let 
Dis-Pat denote the collection of all possible 
distribution patterns, this problem is formulated as 
follows: 

min TCR=^R P X P 

p&DP 

f y Y t p X p =1 Vi e S 

Q T I L-i 1 

^ 1 • < p<=DP 

[ X p e {0,1} Vp e Dis - Pat 

Where R p is the cost rate of distribution 
pattern p e Dis - Pat , and Y? is the binary matrix 
stating whether distribution pattern p ^ Dis -Pat 
covers customer i e S or not. The binary variable 

X p indicates whether distribution 

pattern p e Dis - Pat is selected or not. 

The set Dis-Pat of possible distribution patterns, i.e. 
the number of columns in problem (MP), is infinitely 
large. Therefore, we need an efficient way of 
producing only the most promising distribution 
patterns in Dis-Pat , such that a near-optimal solution 
for the problem is quickly reached. 

Column generation is an appropriate framework to 
this end. We propose column generation procedure 
consists of the following steps: 

Stepl. The set of distribution patterns Dis-Pat is 
initialized with a separate distribution pattern for each 
customer and the initial matrix Y is the identity 
matrix. 

Step2. The linear programming relaxation of the 
master problem (MP) is solved. 

Step 3. Constraint (5) gives a dual price for each of 
the customers, p i ,\/i^S (Generate the basic 
distribution patterns, each serving a single customer 
and use their cost rates as customer priorities p t . Sort 

the customers according to this priority.). These dual 
prices are used as customer priorities. 

Step4 . If the customer priorities are the same as in the 
previous iteration, they are randomized by 
multiplying them with a uniformly generated random 
number between 0 and 10. If this is the sixth time that 
the customer priorities are randomized, the process 
goes to Step 5. 

step (1) Perform the savings heuristic to construct a 
solution and apply the improvement heuristic to it. 
step (2) Perform the insertion heuristic to construct a 
solution and apply the improvement heuristic to it. 

If new distribution patterns are constructed in steps (1) 



and (2), they are added to the set DIS-PAT and the 
corresponding columns are added to the matrix Y. 
The process then returns to Step 2. 

StepS . . The master problem (MP) is solved to 
integrality X p <e {0,1} to provide the final solution. 

In traditional column generation schemes, this means 
that the LP-optimum is found and the column 
generation is finished. However, we are using 
heuristics for the column generating sub problem, 
such that columns with a negative reduced cost may 
not be found even though they do exist, leading to a 

premature ending of the column generation process. 
Another difference with traditional column 
generation is that we add multiple columns per 
iteration. In fact, all distribution patterns that the 
heuristics in the column generating phase encounter, 
are added as new columns in the master problem. For 
large problem instances however, distribution 
patterns are only added to the master problem if their 
reduced cost rate is negative. This is done in order to 
limit the number of columns and thus the CPU-time 
in the Step 5, where a branch-and-bound procedure is 
needed to solve the master problem to integrality. 

4. Comparison to another heuristic 

This section illustrates the value of our solution 
approach which is compared to the solution approach 
of Viswanathan and Mathur [10]. 

As in most papers found in the literature on cyclic 
inventory routing, a fixed vehicle cost is not 
considered in the cost structure of Viswanathan and 
Mathur. They developed a heuristic that generates a 
so-called stationary nested joint replenishment policy 
(SNJRP). They call a policy stationary if 
replenishments are made at equally spaced points in 
time. A nested policy means that if the replenishment 
interval of a given customer is larger than that of 
another customer served by the same vehicle, the 
former is a multiple of the latter. 

Thus, the nestedness corresponds to our delivery 
frequencies, because not all customers visited by the 
same vehicle have the same replenishment interval. 

As a result, the assignment of tours to vehicles is not 
taken into account and the assumption is made that 
there is a separate vehicle for every tour. For the 
cycle times, they only use vehicle capacity as a 
constraining element. 

Compared to our model, the restrictions following 
from (1) traveling, loading and unloading times, and 
(2) limited customer storage capacities are discarded. 
A minimal cycle time for the tours is thus not 
considered. To compare our heuristic method to that 
of Viswanathan and Mathur, the fixed vehicle cost 
has to be discarded and the minimal cycle time T min 
ignored. 

To evaluate the solution approach presented in 
Section3, a number of tests were run on problem 
instances with varying characteristics. The test 
instances are generated according to a 10 x2 3 
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Factorial Design, in which the following three factors 
are considered. 

The first factor is the customer storage capacity 
restriction (CCAP), with level ‘No’ or ‘Yes’ 
indicating the inclusion of the restriction or not. The 
second factor is the holding cost, which is either 80 or 
8 $ per unit per day. This holding cost rate is assumed 
to be the same for all customers. The third and last 
factor is the number of customers (NR). The two 
levels of this factor are [50-90] and [100-150]. This 
factor is a randomized factor, meaning that instead of 
considering two fixed values for the factor, values are 
generated randomly in two distinct intervals. The 
number of customers is thus either between 30 and 70 
or between 80 and 120. 

Comparing our heuristic (Dis-Pat ) to SNJRH as 
Table 1. 



Table 1 Comparing our heuristic (Dis-Pat ) to SNJRH 



HC 


NR 


SNJRH 


Dis-Pat 


A 


100c 


[50-90] 


1420.7 


1289.8 


-9.2 


100c 


[100-150] 


3089.3 


2842.4 


-8.0 


10c 


[50-90] 


721.4 


694.4 


-3.7 


10c 


[100-150] 


1592.2 


1538.1 


-3.4 


Average 




1705.9 


1584.2 


-6.1 



The results shows g that our solution method using 
distribution patterns (Dis-Pat ) and those obtained 
with the stationary nested joint replenishment 
heuristic (SNJRH). Our tests show that this approach 
is more effective than directly solving the problem 
heuristically. 

For all problem instances, solutions are found that are 
cheaper than those proposed by Viswanathan and 
Mathur. On average, there is a cost decrease of no 
less than 6.1%. This is mostly due to the more 
general routing concept that we are using. These 
results show that allowing vehicles to make multiple 
tours gives the opportunity to use vehicle capacity 
much more efficiently and find much better cost 
trade-offs, even without considering fixed vehicle 
costs. Computational experiments show that our 
column generation based heuristic solution approach 
is well capable of finding the appropriate cost 
trade-off under varying circumstances. 

5.Conclusions 

In this paper, we present a practical modeling and 
solution approach for the cyclic inventory routing 
under constant customer demand rates. The proposed 
long-term approach is the first to find three-way cost 
trade-offs between vehicle fleet size costs and 
distribution costs and inventory costs. Furthermore, 
when fleet sizing is not considered, the heuristic 
outperforms an existing heuristic in finding the two 
way trade-off between distribution and inventory 
costs. We assumes a given vehicle type in our 
solution approach. Yet, it can still be applied in cases 
where the fleet contains different vehicle types (with 



different capacities and costs). Actually, if the fleet is 
made of vehicles with different capacities and 
operating costs, then the fourth step in the column 
generation framework must be re-run for each vehicle 
capacity. The mater problem will then select the best 
distribution pattern and thus the appropriate vehicle. 
This of course will increase the computational burden 
of the approach. Also, some other simplifying 
assumptions, including not allowing split delivery 
and imposing a constant time between consecutive 
deliveries may be relaxed, offering more freedom to 
obtain possibly even better solutions. However, this 
will probably increase the complexity of the resulting 
solutions which may, in turn, decrease their 
applicability in practice. 

The cyclic approach brings stability and predictability 
to the planning. On the one hand, this may reduce the 
nervousness and thus the customer demand 
variability and the resulting safety stocks. Thus, 
quantifying this variability reduction is an interesting 
avenue for further research. On the other hand, 
demand variability will not disappear completely, so 
it is worthwhile to investigate how the approach can 
be extended to explicitly take some demand 
variability into account. 
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Abstract 

This paper investigates the use of fuzzy method as a tool 
for model identification and control of a nonlinear and 
multivariable system (greenhouse) when measurement 
data is available. The simulation of the internal climate is 
achieved by a linear model whose coefficients are 
obtained using the fuzzy clustering method. The use of 
fuzzy logic decentralized controller for the regulation of 
climate variables represents powerful way to economise 
the energy and reduce the complexity of the controller. 
Simulation results show that the proposed strategy is 
applied successfully to control the internal climate of 
greenhouse system. 

Keywords: Decentralized Systems, Fuzzy control, 

Greenhouse. 

1. Introduction 

The control of greenhouse climate, in order to improve 
the development of a specific cultivation and to minimize 
the production costs, is becoming increasingly important 
for the growers [9], [12]. The production under 

greenhouse can contribute efficiently to increase the 
productively. To cover these needs it is therefore 
necessary to perfect the air-conditioning of the 
greenhouse in order to maintain the culture in the 
conditions that are compatible with the agriculturist’s 
agronomic and the economic objective [7], [5]. Seen the 
importance of such a process, systems of traditional air 
conditioning used in the habitat (refrigerated machine) 
are too expensive and cannot be put in work in the 
conditions of production. Other methods such as the 
statistical ventilation (roofing), the screens of shadiness 
or the cooling evaporative (moistening) can be adopted. 
Interactions between the internal and external variables 
and complexity of the phenomena (multivariable, 
nonlinear, no stationary) are such that it is often difficult 
to implement these methods. 

Recently, match effort has been devoted to developing 
new control strategies of a greenhouse climate. In this 
respect, a large effort was devoted to develop adequate 



greenhouse climate and crop models, for simulation, 
control and management purposes [6], [11]. The study 
and design of efficient greenhouse environmental 
controllers require having a priori knowledge of the 
greenhouse climate models [10]. These models must be 
related with the external nuances of the outside weather 
conditions (such as solar radiation, outside air 
temperature, wind velocity, etc.), and with the actuating 
actions performed (such as ventilation, cooling, heating, 
among others). The interdependence of the temperature 
and the humidity requires a control strategy which takes 
into account the relationship between these two 
parameters, thus the approach proposed in this work is 
oriented in the synthesis of intelligent climate 
decentralized controller based on the fuzzy logic. 

Fuzzy set theory allows the use of linguistic concepts for 
representing quantitative values and can be employed to 
describe the greenhouse climate based on the system 
identification approach [2], [8]. Compared to traditional 
mathematical modelling, fuzzy modelling possesses some 
distinctive advantage, such as the mechanism of 
reasoning in human understandable terms, the capacity of 
taking linguistic information from human experts and 
combining it with numerical data and the ability of 
approximating complex nonlinear functions with simple 
models. Moreover, fuzzy system modelling is a powerful 
technique to describe complex dynamic systems [4], [5]. 
Based on the results of identification obtained in [1] we 
present, in this work, a fuzzy decentralized controller to 
control the greenhouse air temperature and humidity 
which is a complex system multivariable, nonlinear and 
non stationary. To reduce this complexity we propose to 
decompose the global system MIMO (Multi Input-Multi 
Output) on two subsystems MISO (Multi Input-Single 
Output). 

The paper is organized as follows: in Section 2, we 
describe a greenhouse model. The synthesis of the fuzzy 
decentralized controller is designed in Section 3. In 
Section 4, simulation results are provided to demonstrate 
the effectiveness of the proposed method. 
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2. Formulation problem 

The agriculture greenhouse is an area that permits to 
provide and maintain the creation of a favourable 
environment to plant to growth. The regulation of the 
greenhouse must take in consideration the evolution of 
the meteorological variables and the thermal state. We 
are in presence of MIMO system, non-linear and non- 
stationary in which intervene the energizing exchanges of 
the biologic functions assuring the development of the 
plants. Many works have been done on the development 
of the models of the greenhouse [11], [10], [3], [14] and 
[13] which are grouped on two categories: statistical and 
dynamics models. 

The important number of variables (references, 
perturbations and commands) and the complexity of 
phenomena (biologic, weather, evolution of the plants ...) 
make that we have a multivariable system, non-linear and 
non-stationary [5], (figure l.a). For these reasons, we 
thus preferred to apply a multimodal based on fuzzy logic 
and taken into account all the variables presented on 
figure l.a. 




Figure l.a. Block diagram of the greenhouse model 
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Figure l.b. Input/Output diagram of greenhouse 



where, 7] and H i are respectively temperature and 

relative humidity of the internal air, the perturbation 
variables are T e (external temperature), H e (external 

humidity), R g (solar radiation), V v (wind velocity) and 

the input variable are C h (heating), B r (moistening) and 



O v (roofing). 

To simplify the complexity of the considered system 
greenhouse, we have proposed to decompose the MIMO 
(Multi Input-Multi Output) system into a MISO (Multi 



Input-Single Output) system. By this decomposition, the 
identification and control are collected at level on each 
MISO subsystems; which avoid the computational load 
of the global approach. The nonlinearity, interactions 
between the loops and disturbances are taken into 
account by using the Takagi-Sugeno (TS) modeling. This 
technique can easily transform the nonlinearity problem 
to a linear one. The TS model is defined by the rules, 
which are linear conclusions of the systems outputs 
(internal temperature F and internal hygrometry//.) 

Ti(k + 1) = a l Ov(k) + a 2 Ch{k) + a 3 Br(k) + a 4 Te(k ) 

+ a 5 He(k) + a 6 Rg(k ) + a 7 Vv(k ) + a s Ti(k) + a 9 Hi{k) 

Hi(k + 1) = b x Ov(k) + b 2 Ch(k) + b 3 Br(k) + b 4 Te(k ) 

+ b 5 He(k) + b 6 Rg(k ) + b 7 Vv(k) + b s Ti(k ) + b 9 Hi(k ) 

For the rule j, we obtain: 

If (Ov(k), Ch(k ), . . . , Rd(k)) is Vj then 

Tij (k + 1) = a jX Ov{k) + a j2 Ch(k) + a j3 Br(k ) 

+ a j4 Te(k ) + a j5 He(k ) + a j6 Rg(k ) 

+ a j7 Vv(k) + a jS Ti(k) + a j9 Hi(k) ^2) 



If (Ov(k), Ch(k ), . . . , Rd{k)) is Wj then 
Hij (k + 1) = b Jl Ov(k) + b j2 Ch(k ) + b j3 Br(k) 

+ b j4 Te(k) + b j5 He{k) + b j6 Rg(k) 

+ b j7 Vv(k) + b j% Ti(k) + b j9 Hi(k) 

where 

4 , =K «2 j -■ a 9j] and kj= Uhj K ••• b 9J ] are 
respectively model parameters for temperature and 
humidity for the rule j and 1 < j < M (M is the number 
of rules). 

We define: 

Z(k) = [Ov{k) Ch(k) Br(k) Te(k)He(k) Rg(k) Vv(k) Ti(k) Hi(k )] 

( 3 ) 

and 





, 0, = 


U lj a 2 j 


a 9 j 




J 


K K •• 


• V 



Yj(k) = 0 T j (k)Z(k) 



The vectors y and w. represent a multi dimension 
centers of membership functions [1]. 



3. Fuzzy decentralized control 

The structure of the MIMO fuzzy controller has three 
variables in inputs and two variables in outputs. The 
construction of fuzzy controller is a complex task, 
because many parameters are required for its design. To 
reduce the complexity, we have decomposed the global 
system on two subsystems. The first is defined with two 
inputs and one output, the second is defined with one 
input and one output (figure. 2. a). The decomposition of 
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the system into two subsystems is based on the fact that 
the temperature is much influenced by the heating and 
the action of roofing, the humidity is influenced by the 
moistening. 

For the first subsystem, we define a fuzzy controller, 
named fuzzy controller. 1 (FC1) with four inputs and two 
outputs, for the second subsystem a fuzzy controller.2 
(FC2) with four inputs and one output (figure.2.b). 

Fuzzy controller in closed loop for the greenhouse is 
illustrated in figure. 2. c 

Disturbances 



We use the temperature and humidity variation compared 
to their references eT\ and &H i respectively and their 
variations. 

8 Ti = Tic - Ti , A zTi{k) - zTi(k) - 8 Ti(k - 1) , 
zTe = Te - Ti , A zTe{k) = 8 Te{k) - zTe{k - 1) 
sH / = Hie - Hi , Aeffi(Z) = eHz(£) - sH i(k - 1) , 
zHe = He- Hi , AsH e{k) = sH e(k) - &He(k - 1) . 

All the inputs of the both considered controller FC1 and 
FC2 are three membership functions negative (n), null (z) 
and positive (p) (figure 3): 



Ch k ^ 

OVk 

Commands 



Br k „ 

Command 



Te k He k Rg k Vv k 

Ull 



Sub system 1 

r 



Internal state 
of the temperature 

Ti k 



Interactions 



Sub system 2 



t it t 

Te k He k Ra,Vv k 



_ Hi k 

Internal state 
of the hygrometry 



Disturbances 

Figure 2. a. Inputs and Outputs selected 




with x g {sTi, sTe, sATi, sATe, sHi, sHe, sAZZi, sAZZe} , 
and Lnx, Lzx et Lpx represent the centre of the fuzzy 

subsets A l x , A 2 x et A x . 



Weather disturbances 



STi 

STe 

Asti 

AsTe 



Te k |- 



Rg k Vv k 

. I l 




Ch 

Ov 



Fuzzy controller 1: 

Rj : if £*Ti is A^ Ti and AeTi is A^ AsTi and eTe is A^ Te 
and AsTe is A^ AeTe then Ch is C 1 and Ov is O 1 

( 5 ) 

with 

A k lx g [n,z,p}(k = 1,2,3), C l g {No Heating, Heating), 
O l g [open, no action , closed) and / index of rule. 



Weather disturbances 

Te k He k Rg k Vv k 

l l l l 




Figure 2.b. Fuzzy controllers (FC1 and FC2) 



Disturbances 
Te, He, Rg, Vv 
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He and Tc 



Ti and Hi 



(T\ ^ 


Fuzzy 




System 




Controller 
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Figure 2.c. Fuzzy controller in closed loop for a Greenhouse 



Fuzzy controller 2: 

R) : if sH\ is A^ Hi and Ac Hi is A\ AsHi and sHe is A) sHc 
and AsHe is A\ AeHe then Br is B 1 

( 6 ) 

with B L g {AoMoisterning,Moisterning} 

The goal searched of this work is the regulation of the 
temperature and the humidity inside the greenhouse. In 
the simulation we choose the centre of the fuzzy subsets 
in the order as given in figure 3 :L p =(5,8,3,5,10,10,6,8), 

L z = zeros( 1,8) and L n - —L p 



4 . Simulation results 

In order to validate the proposed approach in simulation 
of day, we have used the results obtained by the classical 
command (on-off) in the presence of the disturbances. 
The file which containing the results of state the day is 
composed of nine colons Ov , Ch , Br , Te , He, Rg, Vv, Ti , Hi 
Using the classical command (on-off), the evolutions of 
the internal, external temperature and hygrometry are 
respectively on figures 4 and 5. 
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Sample 

Figure 4. Temperature curve for the command (on-off) 




Figure 6. Temperature curve for the fuzzy decentralized 
controller 1 




Figure 5. Hygrometry curve for the command (on-off) 



To prove the performance of the proposed approach, we 
have chosen different forms of references. 

First case 

We take the internal temperature and hygrometry 
stemming from the classical command (on-off) as a 
references signal. The figures 6 and 7 show that the 
internal temperature and humidity stemming from the 
proposed fuzzy decentralized controller react correctly 
with the variation of the reference. 

Second case 

The reference signal of temperature and hygrometry are 
fixed (T re f = 20°c, H re f = 80%) for all day. From the 
figures 8 and 9 we can see that a good tracking is 
obtained, which prove the efficiently of the proposed 
approach. 

On the figures 6, 8 and 7, 9 respectively, of the samples 
Oh to 12h and 17h to 24h, we note that the internal 
temperature and hygrometry are closed to the references 
signals, but at 12h to 17h the figures show that there exist 
the oscillations. In the general, the evolution of internal 
temperature and hygrometry are the same at the evolution 
of external temperature and hygrometry with a simple 
amplification. 

Remark: All the oscillations observed on period between 
12h to 17h are naturals, because this period is a critical 
one on the day. 




Sample 

Figure 7. Hygrometry curve for the fuzzy decentralized 
controller 2 




controller 1 (Tref=20°c) . 
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Figure 9. Hygrometry curve for the fuzzy decentralized 
controller l(Href=80%) 



The roofing is active when the temperature is high, this 
corresponds to the activation of the roofing command 
stemming from the classical command (on off). 

We note that the moistening command is not active, the 
same case is recorded for the moistening command 
stemming from the on off command. It was also noted 
that the heating command has not been activated at the 
same time as the classical command (on off), but it was 
active in areas that are characterized by fluctuations 

5. Conclusion 

To reduce the complexity of the greenhouse, witch is 
multivariable, non-linear and non stationary in nature, we 
have proposed to decompose the MIMO system into 
MISO subsystems. By this decomposition, the MISO 
subsystems are represented by compact Takagi-Sugeno 
(TS) model and convert a nonlinear problem into linear 
one. The results obtained, show that the proposed fuzzy 
decentralized controller can be applied to control the 
greenhouse climate. 
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Abstract 

The present paper studies the effect of equipping static 
VAR compensator, SVC in multi-machine power system 
to improve the voltage stability and critical clearing time 
CCT. of multimachine power system. Three static VAR 
compensators, SVCs are connected with certain power 
system buses to stabilize the system voltage and improve 
the overall power system stability subjected to different 
disturbances for a wide range of operating conditions . 
Also, the SVCs improve the power system critical 
clearing time to a significant values in case of unrepeated 
three phase short circuit fault at all system buses for 
different operating conditions. The results indicate the 
importance of installing SVCs to grid power system due 
to their effect not only on stabilizing the system voltages 
but also increasing the power system critical clearing 
time. The critical clearing time is obtained by increasing 
the fault time interval until the system loses its stability. 
Increasing the CCT can help in the reliability of the 
protection system and reduce the protection system rating 
and its cost. The studied 3 -machine nine- bus power 
system is modeled and simulated using MATLAB 
software package. 

Keywords 

Multi-machine power system - Voltage stability 
improvement - Static VAR Compensator - Critical 
clearing time. 



Nomenclature 

5i(t) Generator rotor angle 



Relative speed of the i gen. in p.u. 
Synchronous speed in p.u. 

Mech. input torque, in p.u.. 

Elec, torque, in p.u. . 

Transient EMF in quad, axis, in p.u.. 
Direct axis current, in p.u 
Quad, axis current, in p.u... 

Per unit damping constant. 

Inertia constant in second. 



(Di 

(D s 
T 

1 m 

T e i 
Eq 

Idi 

Iqi 

Di 
Hi 

T doi Direct axis transient s. c. time const, in sec. 



Tqoi Quad, axis transient s. c. time const, in sec. 

x doi Direct axis reactance in p.u. 

x do i Quad, axis reactance in p.u. 

x doi ^ Direct axis transient reactance in p.u. 

x qoi v Quad, axis transient reactance in p.u. 

Eg Field voltage in per unit. 

K Ai Gain of the excitation amp. 

B l : Inductive susceptnce of the SVC. 

B c : Capacitive susceptnce of the SVC 

1. Introduction 

Flexible Alternating Current Transmission System 
(FACTS) devices are used extensively to increase the 
amount of power transferred of the transmission systems 
and increase the power system stabilization 1-4]. Static 
VAR compensators (SVCs) are used in power systems in 
both transmissions and distribution systems. The main 
objective of installing SVCs in electric power system is 
to provide rapid control of voltages at weak points in 
electric power systems. SVCs is used also to provide 
system stabilization in case of insufficient damped inter- 
area oscillations and improve the power system power 
quality [5-6]. SVCs should be controlled in order to 
provide rapid and continuous reactive power during 
steady state and dynamic operating conditions. The 
application of SVCs in Extra High Voltage, EHV 
network to provide a dynamic reactive power 
compensation to keep the bus voltages within the 
permissible values under varying load conditions, and 
improve voltage stability. The location and size of SVCs 
have been discussed in details in both transient and 
steady state [7]. The application of multiple compensation 
have been considered by the utilities to meet the 
performance objectives. These objectives are considered 
based on evaluate different goals, such as improving 
power system network voltage profile, reducing the 
power system losses, increasing the power transmissions, 
damping the power system oscillations, and improving 
the overall power system stability [8- 13]. Recently, there 
are different control techniques have been applied to 
provide (SVCs) rapid response such as self tuning, fuzzy 
logic, neural network, neurofuzzy[14-17]. 
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Figure 1 IEEE WSCC 3 Machine 9-bus Power System 



The stability analysis of power systems may consider the 
determination of the system critical clearing time, CCT 
for a given fault in order to find the maximum value of 
CCT for which the system is still stable. If the fault is 
cleared within this time the system will remain stable. 
However if the fault cleared after this time the power 
system will lose its stability. The calculation of CCT is 
very important from the protection point of view. This 
will aid to determine the characteristics of protection 
system required to the power systems. There are several 
conventional methods to compute the CCT, such as 
numerical integration and Lyapunove- 
type stability criteria[18]. Recently the neural network 
approach have been used extensively to determine the 
critical clearing time of power systems [19-23]. 

This paper presents the effect of installing static VAR 
compensators (SVCs) to multi-machine power systems. 
The results obtained investigate that installing SVCs 
stabilize the power systems voltages and damp the 
oscillations in case of disturbance very fast. The studied 
multi-machine power system have been simulated and 
the critical clearing time is obtained with and without 
installing (SVCs). The result indicates that with installing 
(SVCs) the critical clearing time, CCT have been 
increased. The studied system is IEEE WSCC 3 
Machines 9-bus Power System. 



2. Modeling of the Presented Power System 

The multimachine power system under study with three 
SVCs connected at certain buses is shown in Figure 1. 
Different models have been used for stability studies. 7 th 
order model for the time simulation of studied power 
system is implemented . Each synchronous generator is 
equipped with IEEE -Type I Exciter model[24]. 



= CO; -<B q 

dt 1 s 


(1) 
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-(Ela+X^-D^-rtO 
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(5) 
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The algebraic equations are: 

V ti = ( V di 2 + V qi 2 ) m 


(8) 


Vdi =E d i - Re Idi + X e I q i 


(9) 


Vqi — Eqi - R e Iqi — X e I^i 


(10) 


I=YV 


(11) 


Ii=(id+ji q )e j(5 i -’ t/2) 


(12) 


V,=(V d +jv q ) e i(8 i ' lt/2) 


(13) 


V. ;i — Vj e|5i- 


(14) 



For the study of the system stability, the previous 
equations are solved by Matlab software. 

3. Static VAR Compensator Model 

The SVC model used in this study is fixed capacitor, 
thyristor controlled reactor(FC-TCR) type as shown in 
figure 2. By changing the thyristor firing angle 
(alpha), the reactor current can be controlled according to 
the following equations [14 ] : 




The SVC model is added to the system equations and 
solved using matlab software to obtain the critical 
clearing time CCT of the studied power system and study 
the effect of installing SVC in voltage stability 
enhancement. 

4. PI Controller 

PI controller is used in this study due to it's applicable 
device. The block diagram of the SVC PI controller is 
shown in figure 3. The equation describe the PI controller 
can be expressed as: 

a (t) = AT [Aco(t)Ki + K P (Aoo(t) - A co (t-1)))] + a (t-1) 

( 20 ) 

Where: ti/ 2 < a (t) < n 

H max 




Figure 2 Schematic Diagram of FC-TCR SVC type 



5. Results and Discussions 

The parameter data of studied power system are given in 
the appendices . The load flow study is performed firstly 
to determine the critical clearing time of the studied 
power system. Different load conditions are applied in 
order to evaluate the effect of installing SVC on the 
power system critical clearing time and the bus voltages 
after applying mechanical input power disturbance. 
Table 1 gives the steady state load flow results of the 
studied power at different operating conditions. Figure 7 
shows the flowchart of computer solution to the studied 
system. 



h(a) = 



V (In - 2a + sin 2 a) 
ncoL 



where < a < 71 
cr- sin cr 

h(v) = 7 — V 

ncoL 

The reactor susceptance is given as : 

h ( ct ) 



B L (cr)-- 



V 



(15) 



(16) 



(17) 



From equations (16),(17) the susceptance can be obtained 
as : 

cr - sin cr 



b l O) = ■ 



ncoL 



(18) 



However the equivalent SVC susceptance is given by the 
summation of the fixed capacitor susceptance and the 
thyristor controlled reactor susceptance as follows: 

^ svc ~ — B L (cr) (19) 

The conduction angle a and the firing angle a are related 
by the following equation: 
cr = 2(n - a) 



Table 1 : The buses voltages of the Load flow results of the 
s tudied system at different operating conditions. 



Load % 
Bus # 


25 


50 


75 


100 


125 


150 


1 


1.040 


1.040 


1.040 


1.040 


1.040 


1.040 


2 


1.025 


1.025 


1.025 


1.025 


1.025 


1.025 


3 


1.025 


1.025 


1.025 


1.025 


1.025 


1.025 


4 


1.057 


1.048 


1.037 


1.026 


1.013 


0.998 


5 


1.059 


1.040 


1.019 


0.996 


0.969 


0.939 


6 


1.064 


1.049 


1.032 


1.013 


0.992 


0.968 


7 


1.050 


1.043 


1.035 


1.026 


1.014 


1.001 


8 


1.053 


1.042 


1.030 


1.016 


1.000 


0.981 


9 


1.051 


1.046 


1.039 


1.032 


1.024 


1.014 



5.1 Critical clearing time evaluation 

In order to obtain the critical clearing time for the studied 
power system , the power system load is increased 
gradually and load flow study is executed for the stability 
study. The CCT can be obtained by increasing the time 
interval of the fault gradually till the studied system 
loses its stability. To study the powerful of the proposed 
SVC controller, three SVCs are connected at the buses 
number 4, 7, 9. of the studied power system as shown in 
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figure 1 and the fault interval time increases gradually up 
till the system loses its stability and then the CCT can be 
determined . 

Table 2 to Table 10. shows the studied power system 
critical clearing time with and without installing static 
VAR compensator at different load level. It is very clear 
that the effect of installing SVC in increasing the CCT up 
to 200% in some cases very clear. This will aid to choose 
the suitable protection system. 



Table [2]: Critical clearing time with and without the proposed SVC at 
di fferent loading levels when the Fault occurs at bus # 1 : 



% of load level 


CCT without SVC 
In m.sec. 


CCT with SVC 
In m.sec 


25 


122 


223 


50 


116 


215 


75 


110 


193 


100 


104 


189 


125 


98 


182 


150 


92 


176 



Table [3]: Critical clearing time with and without the proposed SVC at 
differ ent loading levels when the Fault occurs at bus # 2: 



% of load level 


CCT without SVC 
In m.sec. 


CCT with SVC 
In m.sec 


25 


133 


249 


50 


126 


242 


75 


119 


233 


100 


108 


216 


125 


101 


208 


150 


95 


187 



Table [4}: Critical clearing time with and without the proposed SVC at 
differ ent loading levels when the Fault occurs at bus # 3: 



% of load level 


CCT without SVC 
In m.sec. 


CCT with SVC 
In m.sec 


25 


139 


269 


50 


131 


258 


75 


121 


247 


100 


114 


223 


125 


109 


205 


150 


101 


195 



Table[5]: Critical clearing time with and without the proposed SVC at 
differ ent loading levels when the Fault occurs at bus #4: 



:% of load level 


CCT without SVC 
In m.sec. 


CCT with SVC 
In m.sec 


25 


434 


845 


50 


401 


823 


75 


369 


673 


100 


256 


494 


125 


109 


217 


150 


91 


163 



Table [6]: Critical clearing time with and without the proposed SVC at 
differ ent loading levels when the Fault occurs at bus # 5 



% of load level 


CCT without SVC 
In m.sec. 


CCT with SVC 
In m.sec 


25 


161 


307 


50 


126 


231 


75 


105 


186 


100 


91 


147 


125 


82 


123 


150 


73 


97 



Table [7}:Critical clearing time with and without the proposed SVC at 
differ ent loading levels when the Fault occurs at bus #6: 



% of load level 


CCT without SVC 
In m.sec. 


CCT with SVC 
In m.sec 


25 


134 


253 


50 


127 


238 


75 


122 


227 


100 


109 


221 


125 


103 


213 


150 


98 


202 



Table [8]:Critical clearing time with and without the proposed SVC at 
differ ent loading levels when the Fault occurs at bus # 7: 



% of load level 


CCT without SVC 
In m.sec. 


CCT with SVC 
In m.sec 


25 


462 


834 


50 


421 


713 


75 


343 


587 


100 


223 


386 


125 


141 


237 


150 


123 


203 



Table [9]:Critical clearing time with and without the proposed SVC at 
differ ent loading levels when the Fault occurs at bus # 8: 



% of load level 


CCT without SVC 
In m.sec. 


CCT with SVC 
In m.sec 


25 


503 


969 


50 


482 


933 


75 


463 


902 


100 


441 


857 


125 


234 


441 


150 


141 


265 



Table [10]: Critical clearing time with and without the proposed SVC at 
differ ent loading levels when the Fault occurs at bus # 9: 



% of load level 


CCT without SVC 
In m.sec. 


CCT with SVC 
In m.sec 


25 


471 


912 


50 


451 


881 


75 


437 


867 


100 


221 


408 


125 


147 


267 


150 


113 


203 



5.2 Voltage Stability Study 

In order to indicate the effect of installing SVCs with 
each synchronous machine on enhancement of power 
system voltage stability, 0.02 p.u. input mechanical 
power disturbance is considered to each machine 
separately. Figure 4 depicts the power system bus voltage 
response with and without static VAR compensator when 
2% input mechanical power disturbance is applied at the 
three machines each separately when the load is 50%. 
The effect of using SVC is very clear in damping the 
voltage oscillations very fast and improve the system 
performance. Figure 5 shows the machine bus voltages 
when the load is exceed to full load and 0.02 p.u. input 
mechanical power is considered. The performance with 
installing SVC is very clear in damping the bus voltages 
oscillations very fast. Figure 6 depicts bus voltages when 
the studied power system load became 125% with and 
without installing SVC. The effect of SVC in damping 
the machine bus voltages oscillation is clear. 
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Time in sec. 



Figure 4 Studied System Voltage when 0.02 p.u. increased in 
input mechanical power (50% load) 



Figure 6 Studied System Voltage when 0.02 p.u. increase in 
input mechanical power (125% load) 



1.042 
1.041 | 
.E 1.04 


l \ ) 


> 
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With SVC 


1 038q 


i 1 23456789 10 

Time in sec. 


1.027 

1.026 


J A 


Vb2 ii 


j 
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0 


i 1 23456789 10 

Time in sec. 


1.0255 

-i 1 025 


| SE j| f J; 

i 1 V P £ 


Vb3 in p.i 


if 11 M o u f 

I 1 V . Without SVC 

With SVC 



1.0235 1 1 1 1 1 1 1 1 1 

01 23456789 10 



Figure 5 Studied System Voltage when 0.02 p.u. increase in 
input mechanical power (100% load) 



6. Conclusions 

This paper investigates the importance of equipping static 
VAR compensator with the power systems. Different 
load levels and different fault locations of the studied 
power system are suggested to study the effect of 
installing SVC in increasing the CCT, and stabilizing the 
studied power system voltages. The simulation results 
proved that connecting of SVC to the power system 
increases the power system critical clearing time with a 
specific values when a three phase short circuit fault is 
occurred. The installing SVCs to power systems can 
increase the CCT to very close to 200% in compared 
without installing SVCs. Also the power system voltages 
became more stable very fast when 0.02 p.u. of an input 
mechanical disturbance is considered to the synchronous 
machines at different operating conditions. 
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Appendix 

The 3 machine nine bus power system line data 



Line 


SB 


EB 


Line 


Ysh 


Zs 


g 


jb 


r 


X 


1 


1 


4 


1 


0.00 


0.000 


0.000 


0.057 


2 


4 


5 


1 


0.00 


0.088 


0.010 


0.085 


3 


4 


6 


1 


0.00 


0.092 


0.017 


0.092 


4 


5 


7 


1 


0.00 


0.153 


0.032 


0.161 


5 


6 


9 


1 


0.00 


0.179 


0.039 


0.170 


6 


9 


3 


1 


0.00 


0.000 


0.000 


0.058 


7 


7 


2 


1 


0.00 


0.000 


0.000 


0.062 


8 


9 


8 


1 


0.00 


0.104 


0.011 


0.100 


9 


8 


7 


1 


0.00. 


0.074. 


0.008. 


0.072. 



Three machines data: 





Mach. 1 


Mach.2 


Mach. 3 


H 


23.64 


6.4 


3.01 


x d (p.u) 


0.1460 


0.8958 


1.3125 


x d (p.u) 


0.0608 


0.1198 


0.1813 


Xq (p.u) 


0.0969 


0.8645 


1.2578 


x Q ' p.u.) 


0.0969 


0.1969 


0.2500 


Tdo sec. 


8.960 


6.00 


5.89 


T q sec. 


0.310 


0.535 


0.600 
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Figure 7 Flowchart of solution the studied system 



Exciter data: 





K a 


T A (sec.) 


K E 


T E (sec.) 


k f 


T F (sec.) 


Excit. 1 


20 


0.2 


1.0 


0.314 


0.63 


0.35 


Excit. 2 


20 


0.2 


1.0 


0.314 


0.63 


0.35 


Excit. 3 


20 


0.2 


1.0 


0.314 


0.63 


0.35 



Power System Loads, in KVA: 
Load A = 125 + j 40, LoadB = 
Load C = 100 + j 30. 



90 +j 30, 
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Abstract 

Optimum design of Switched Reluctance Motor (SRM) 
requires knowledge of its parameters, such as inductance, 
flux linkage and torque. In this paper, the model of SRM 
with two different stator pole shapes have been 
developed and analyzed with ANSYS (Analysis System) 
Workbench. The variation of inductance, flux density and 
directional force with respect to rotor position for both 
the models have been gathered, which are the 
fundamental characteristics of SRM and have significant 
influences on motor performance. It is observed that the 
proposed parallel slot type stator pole (model 2) gives the 
improved flux density and directional force/ Torque. 

Keyword: Switched Reluctance Motor, ANSYS, Flux 
linkage, Directional Force/Torque. 

1. Introduction 

The theory of motor with variable reluctance has been 
started to be used since 1980’s for the variable or 
adjusted speed application and the use of these motors 
started to be very popular in engineering applications. 
The advantages brought by the cheap and powerful 
switching elements enabled the re-discovery of this 
motor [1]. 

Many advanced manufacturing processes require smooth, 
high precision and/or high speed robotic manipulation in 
3D space. Electrical Motors, as energy conversion 
devices, are widely used in industrial applications. The 
simulation and experimental investigations of the 
electromagnetic field distribution are one of critical topic 
for developing electromagnetic products [2]. Due to 
exclusive features of the SRM such as lack of any coil or 
permanent magnet on the rotor, simple structure and high 
reliability, this makes it a suitable candidate for operation 
in hard or sensitive conditions. Although the FEM was a 
time-consuming and difficult method in the past, with 
contemporary high speed computers with large memory, 
the method can be used properly in analysis of electric 
machines. In addition, various commercial finite-element 
packages ANSYS Workbench have been so far 
developed for this purpose [3]. 



In this paper, Switched Reluctance motor with two 
different pole configurations were taken for magnetic 
analysis. Section (2) describes about the parametric 
modeling of motors. Finite element method and 
assumptions made for the analysis were described in 
Section (3). Finally, a comparison of the two structures 
and conclusions are presented in Section (4) and (5) 
respectively. 

2. The Parametric Simulation Model 

The ANSYS Workbench platform is a powerful multi- 
domain simulation environment. Multiphysics simulation 
from ANSYS provides powerful simulation tools for 
solving industry’s toughest multiphysics challenges. It 
also includes solutions for both direct and sequentially 
coupled physics problems including direct coupled-field 
elements [4]. 

In order to analyze the SRM by ANSYS Workbench, the 
geometrical model of the motor was plotted at first and 
appropriate materials are assigned to different regions [5]. 
The SRM has a complicated configuration and it is 
cumbersome and time-consuming to draw the 
geometrical model every time before simulation. This is 
not practical when the simulation model is supposed to 
be used in optimal design of the SRM where the design 
process should be repeated frequently. Therefore, a 
parametric geometric model is desirable [2,3]. 

With a reasonable approximate, the whole geometric 
structure of SRM is generated by geometric parameters, 
including radius of different parts of the motor, stator and 
rotor pole arc lengths. The design procedure for SRM is 
referred from [6, 7]. The dimensional details of 100 W, 
3200 rpm SRM with 6/4 pole configuration is shown in 
the Table 1. and the cross section is shown in Figure 1. 
The details of stator pole configuration of model 1 is 
shown in Figure 2 and Table 2 respectively. Similarly, 
for the proposed model (model 2) details are shown in 
Figure 3 and Table 3 [7, 8]. 
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Figure 1. 6/4 Switched Reluctance Motor Geometry 



Table 1 : The Main Dimensions of SR Motor 



Symbol 


Meaning and Unit 


Value 


Dso 


Stator outer diameter (mm) 


122.0 


Dsi 


Stator inner diameter (mm) 


106.0 


Dro 


Rotor outer diameter (mm) 


65.4 


Dri 


Rotor inner diameter (mm) 


46.0 


Dsf 


Shaft diameter (mm) 


28.0 


lg 


Air gap (mm) 


0.3 


Lstk 


Stack length (mm) 


15.0 


ns 


Stator pole number 


6 


nr 


Rotor pole number 


4 


q 


Number of phase 


3 




Table 2: Details of Stator Pole Shapes of Model 1 



Symbol 


Meaning and Unit 


Value 


Twd 


Tooth Width (mm) 


13 


so 


Slot Open Length (mm) 


15 


Tag 


Tooth Tang Angle (degree) 


30 




Table 3: Details of Stator Pole Shapes of Model 2 



Symbol 


Meaning and Unit 


Value 


Twd 


Tooth Width (mm) 


13 


SO 


Slot Open Length (mm) 


15 


Tag 


Tooth Tang Angle (degree) 


30 


Swd 


Slot Width (mm) 


28 



3. Finite Element Model 



The magnetic field inside the motor is determined by 
computing the magnetic vector potential A. This 
satisfies the non-linear Poisson’s equation. 



d_ 

dx 



dA } d 



- +- 



dx ) dy 



dA 

dy 



dz 



y- 



aA 

dz 



= -J 



( 1 ) 



in a three-dimensional Cartesian system where y is the 
magnetic reluctivity and J is the current density vector. 
The energy functional is a generalized from of a linear 
technique [8]. 

B _ _ a _ 

F = J J H.dBdv - J J J.dAdv (2) 

VO VO 

where B is magnetic flux density in Tesla (T), H is 
magnetic field intensity in Amps per meter (A/m), J is 
Current density in Amps per square meter (A/m 2 ). 

The 3D Finite Element Method (FEM) Analysis of SR 
motor is carried out with help of ANSYS workbench 
with the following assumptions 

1 . The current density within any conductor is uniform. 

2. Hysteresis effects are neglected. 

3. The induced conduction current in iron is neglected 
because of the relatively high resistance offered by 
steel laminations for eddy currents in the axial 
direction. 

The simulation and experimental investigations of the 
electromagnetic field distribution of SR motor are 
analyzed with help of what-if analysis for different 
rotor positions for the given excitation current. 

Torque characteristic is one of the most important 
electromagnetic characteristics of SRM, which can be 
used in calculation of average torque and torque 
ripple. [9, 10] In fact, this static characteristic shows 
torque exerting on the rotor for different rotor 
positions. Direct calculation of torque in SRM with 
analytical method is impossible. Usually, torque is 
derived from the flux-linkage characteristic of phase 
winding. However, analytical methods cannot model 
precisely the region in where the stator and rotor poles 
begin to overlap, so prediction of the static torque is 
faced with a problem in this region [10,12-14]. 

Using ANSYS Workbench, it is possible to predict the 
torque directly. It is just enough to exert Maxwell work 
torque boundary conditions on the rotor [3]. 



4. Results and Discussion 

• The variation of inductance for various rotor 
positions is shown in Figure.4. For model 2 it is 
found that in aligned position the variation of self 
inductance is minimum. 

• The flux density variation for the both model is 
shown in Figure. 5. In the model 2 variation of flux 
density is found less, will give more torque. 
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• Comparison of Figure. 7 and Figure. 8 gives the 
result of improvement of total force on rotor in 
model 2. 

• The variation of directional force/torque along the 
X axis is shown in Figure. 9 and Figure. 10 it 
indicates that model 2 will give more torque than 
model 1. 

• The Directional force/torque along y-axis from 
Figure. 11 and Figure. 12 can be found more in 
model-2. 

• The flux density variation along motor body for 
both models is shown in Figure. 13 and Figure. 14. 
in model 2 it is clear that ‘B’ is increased. 




25.00 



Figure. 7. Total Force in Rotor for Model 1 










Figure. 4. Self Inductance of winding at various Rotor positions 



Total Force 

Type: Total Force 
Unit: N 
Time: 2 

2/26/2009 12:28 PM 

0.33978 Max 

0.30203 
0.26427 
0.22652 
J 0.18877 
■ 0.15101 
0.11326 
0.075507 
0.037753 
1.4064e-9 Min 



4 



x 





* * * * & & A'V * 

Rotor Angle (Dcsrccs) 



Figure. 5. Flux Density Comparison at various Rotor positions 



Figure. 8. Total Force in Rotor for Model 2 




25.00 



Figure. 9. Directional Force/Torque along X-axis in Rotor for Model 1 




Rotor Angle (Dcarccs) 




30.00 



Figure. 6. Total force Comparison at various rotor positions 



Figure. 10. Directional Force/Torque along X-axis in Rotor for model 2 
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Directional Force Torque 2 

Type: Directional Force/Torque (YAxis ) 
Unit: N 
Time: 2 

2/1 9/2009 8:58 AM 



Figure. 1 1 . Directional Force/Torque along Y-axis in Rotor for model 1 



Directional Force/Torque 2 

Type: Directional Force/Torque (YAxis) 
Unit: N 
Time: 2 

2/1 9/2009 1 0:08 AM 



1.1986 Ma! 

0.92783 
-I 0.65709 

- 0.38635 
m o n 562 
■ -0.15512 

- -0.42586 
-| -0.69659 

-0.96733 
-1.2381 Mir 







Figure. 12. Directional Force/Torque along Y-axis in Rotor for Model 2 




Total Flux Density 

Type: Total Flux Density 
Unit: mT 
Time: 2 

2/1 9/2009 8:59 AM 

1208 Max 

1073.8 

- 939.58 

- 805.35 

- 671.13 
3 536.9 j 

402.61 

- 268.4! 

134.2: 



Figure. 13. Variation of Flux Density for Model 1 




Figure. 14. Variation of Flux Density for Model 2 



5. Conclusion 

In this paper three-dimensional magnetic analysis of a 
Switched Reluctance Motor with two different pole 
configurations have been given. From the analysis the 
variation of inductance, flux density & directional force / 
torque of both the model have been predicted. The results 
shows that the proposed pole configuration (model 2) is 
advantageous than the existing pole configuration 
(model l).This study can be used for further modification 
of physical size of the motor to get more improved 
performance and to manufacture the SRM. In future the 
optimum shape and dimensions of SRM stator and rotor 
pole can be identified using soft computing techniques. 
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